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Context. Emission lines from protoplanetary disks originate mainly from the irradiated surface layers, where the gas is generally 

■ warmer than the dust. Therefore, the interpretation of emission lines requires detailed thermo-chemical models, which are essential to 
convert line observations into understanding disk physics. 

Aims. We aim at hydrostatic disk models that are valid from 0.1 AU to 1000 AU to interpret gas emission lines from UV to sub-mm. 
, In particular, our interest lies in the interpretation of far IR gas emission lines as will be observed by the Herschel satellite, related to 

^ I' the Gasps open time key program. This paper introduces a new disk code called ProDiMo. 

Methods. We combine frequency-dependent 2D dust continuum radiative transfer, kinetic gas-phase and UV photo-chemistry, ice 
formation, and detailed non-LTE heating & cooling with the consistent calculation of the hydrostatic disk structure. We include 
Fell and CO ro- vibrational line heating/cooling relevant for the high-density gas close to the star, and apply a modified escape 
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ABSTRACT 



o 

^ . probability treatment. The models are characterized by a high degree of consistency between the various physical, chemical and 

^ I ' radiative processes, where the mutual feedbacks are solved iteratively. 

Results. In application to a T Tauri disk extending from 0.5 AU to 500 AU, the models show that the dense, shielded and cold midplane 
(z/rS 0. 1, Tg » Tii) is surrounded by a layer of hot (Tg » 5000 K) and thin («<h> ~ 10' to 10* cm"^^) atomic gas which extends radially 
^ ■ to about 10 AU, and vertically up to z/r x 0.5. This layer is predominantly heated by the stellar UV (e. g. PAH-heating) and cools 

' via Fell semi-forbidden and Oi630nm optical line emission. The dust grains in this "halo" scatter the star light back onto the disk 

' which impacts the photo-chemistry. The more distant regions are characterized by a cooler flaring structure. Beyond 100 AU, Tg 

^■f^ , decouples from even in the midplane and reaches values of about Tg a 27^. 

. Conclusions. Our models show that the gas energy balance is the key to understand the vertical disk structure. Models calculated 

with the assumption Tg = show a much flatter disk structure. The conditions in the close regions (< 10 AU) with densities n^^^y « 10* 
to 10'^ cm"' resemble those of cool stellar atmospheres and, thus, the heating and cooling is more stellar-atmosphere-like. The 
application of heating and cooling rates known from PDR and interstellar cloud research alone can be misleading here and more work 
needs to be invested to identify the leading heating and cooling processes. 
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?H ' 1 ■ Introduction but still provides enough ionization to drive a rich ion-molecule 

F~ chemistry (Bergin et al., 2007). However, the sensitivity of cur- 

The structure and composition of protoplanetary disks play a key ^ent radio telescopes allows only observations of small samples 

role in understanding the process of planet formation. From ther- (^gnt 2005) and in a few cases detailed studies of individ- 

mal and scattered hght observations, we know that protoplane- ^^1 objects (e.g. Qi et al., 2003; Semenov et al., 2005; Qi et al., 

tary disks are ubiquitous in star forming regions and that the 2OO8). The gas temperature in the disk surface down to contin- 

dust in these disks evolves on timescales of 10" yr (Haisch et al., ^^j^ optical depth of ~ 1 decouples from the dust temperature 

2006; Watson et al., 2007). However, dust grains represent only ^nd ranges from a few thousand to a few hundred Kelvin (Kamp 

about 1% of the mass m these disks - 99% of their mass is gas. & Dullemond, 2004; Jonkheid et al., 2004; Dullemond et al.. 

While observations of the dust in these systems have a long 2007). The hot inner disk can show ro-vibrational line emission 

history (e.g. Beckwith et al., 1990), gas observations are intrin- in the near-IR either due to fluorescence (at the disk surface) 

sically more difficult and have focused until recently on rota- and/or thermal excitation (e.g. Bary et al., 2003; Brittain et al., 

tional lines of abundant molecules such as CO, HCN, HCO+ 2007; Bitner et al., 2007). More recently, near-IR gas lines have 

(e.g. Beckwith et al., 1986; Koerner et al., 1993; Dutrey et al., also been detected in Spitzer IRS spectra, revealing the pres- 

1997; van Zadelhoff et al., 2001; Thi et al., 2004). These lines ence of water, H2 and the importance of X-rays (Pascucci et al., 

generally ti-ace the outer cooler regions of protoplanetary disks 2007; Lahuis et al., 2007; Salyk et al., 2008). The launch of the 

and probe a layer at intermediate heights, where the stellar- UV Herschel satellite in 2009, opens yet another window to study 

radiation is sufficiently shielded to suppress photodissociation, the gas component of protoplanetary disks through the domi- 

nant cooling lines [Oi], [Cii] at the disk surface as well as many 

* The Scottish Universities Physics Alliance additional molecular tracers of the warmer inner disk such as 
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water and CO. The study of the gas in protoplanetary disks is 
the main topic of the Herschel open time Key Program "Gas 
in Protoplanetary Systems" (Gasps, PI: Dent). Other guaranteed 
and open time Key Programs, such as "Water in Star Forming 
Regions with Herschel" (Wish, PI: van Dishoeck) and "HIFI 
Spectral Surveys of Star Forming Regions" (PI: Ceccarelli), will 
also observe gas lines in a few disks. 

Disk structure modeling was initially driven by dust obser- 
vations and developed from a simple two-layer disk model (e.g. 
Chiang & Goldreich, 1997) into detailed dust continuum ra- 
diative transfer models that are coupled with hydrostatic equi- 
librium (e.g. D'Alessio et al., 1998; Dullemond et al., 2002; 
Dullemond & Dominik, 2004; Pinte et al., 2006). The assump- 
tion in all these models is that gas and dust are well coupled and 
the hydrostatic scale height then follows from the dust tempera- 
ture. However, the gas temperature decouples from the dust tem- 
perature and the vertical disk structure will adjust to the gas scale 
height, forcing the dust to follow if it is dynamically coupled. 
This approach has been followed by Nomura & Millar (2005) 
and Gorti & Hollenbach (2004, 2008). Nomura & Millar use a 
small chemical network (only CO, C^ and O) and a limited num- 
ber of heating/cooling processes namely photoelectric heating, 
gas-grain collisions and line cooling from [Oi], [Cii] and CO. 
Gorti & Hollenbach use an extended set of reactions (84 species, 
~ 600 reactions) and the relevant low-density heating/cooling 
processes drawn from photo dissociation region (PDR) physics. 
Other models do not solve for the vertical hydrostatic disk struc- 
ture (e. g. Kamp & Dullemond, 2004; Meijerink et al., 2008; 
Woods & Willacy, 2008). Kamp & Dullemond use a chemical 
reaction network of ~ 250 reactions among 48 species and a set 
of heating/cooling processes comparable to Gorti & Hollenbach 
(2004). The models of Meijerink et al. focus entirely on the X- 
ray irradiation of the disk, thus excluding UV processes; the 
chemical reaction network is limited to 25 species and 125 re- 
actions. Woods & Willacy use again a standard set of PDR heat- 
ing/cooling processes, but also account for X-rays. Their chemi- 
cal network includes 475 gas and ice species connected through 
~ 8000 gas phase and surface reactions. 

This paper presents a new disk code that includes additional 
heating/cooling processes relevant for the high densities and 
high temperatures present in the inner parts of the disk, resem- 
bling the conditions in tenuous atmospheres of cool stars. The 
models are characterized by a high degree of consistency be- 
tween the various physical, chemical and radiative processes. In 
particular, the results of a full 2D dust continuum radiative trans- 
fer are used as input for the UV photo-processes and as radiation 
background for the non-LTE modelling of atoms and molecules 
to calculate the line heating and cooling rates. This allows the 
models to extend closer to the star and include modelling of the 
so-called inner rim. 

The paper is structured as follows. Section 2 introduces the 
new code ProDiMo and presents the concept of global iterations. 
Section 3 describes the assumptions used to calculate the hy- 
drostatic disk structure including "soft edges". In Sect. 4, we 
present the 2D dust continuum radiative transfer with scattering 
and band-mean opacities. Section 5 summarizes the gas-phase 
and photo-chemistry dependent on the UV continuum transfer 
results. In Sect. 6, we outline the heating and cooling rates in- 
cluded in our model and present a modified escape probability 
method. Section 7 closes the theory part of the paper with the 
calculation of the sound speeds as preparation of the next calcu- 
lation the the disk structure. We apply ProDiMo to a standard 
TTauri-type protoplanetary disk with disk mass 0.0 IM© which 
extends from 0.5 AU to 500 AU in Sect. 8. The resulting phys- 
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Fig. 1. Concept of global iterations in ProDiMo. The circular arrows 
on the r.h.s. indicate sub-iterations. For example, the dust temperature 
structure needs to be iterated in the continuum radiative transfer. 



ical and chemical structure of the disk is shown and compared 
to a model where we assume T„ - T^. We conclude the paper in 
Sect. 9 with an outlook to future applications. 

2. ProDiMo 

ProDiMo is an acronym for Pro toplanetary Disk Model. It is 
based on the thermo-chemical models of Inga Kamp (Kamp 
& Bertoldi, 2000; Kamp & van Zadelhoff, 2001; Kamp & 
Dullemond, 2004), but completely re-written to be more flexi- 
ble and to include more physical processes. 

ProDiMo uses global iterations to consistently calculate the 
physical, thermal and chemical structure of protoplanetary disks. 
The iterations involve 2D dust continuum radiative transfer, gas- 
phase and photo-chemistry, thermal energy balance of the gas, 
and the calculation of the hydrostatic disk structure in axial sym- 
metry (see Fig. 1). The different components will be explained 
separately in the forthcoming sections. 

Physical processes not yet included are X-ray heating, X- 
ray chemistry, spatially dependent dust properties, and PAH- 
chemistry. These processes will be addressed in future pa- 
pers. ProDiMo is under current development. The code can 
be downloaded from https://forge.roe.ac.uk/trac/ProDiMo, start 
at https://forge.roe.ac.uk/trac/ROEforge/wiki/NewUserForm to 
get a ProDiMo user account. 



3. Hydrostatic disk structure 

We consider the hydrostatic equation of motion in axial symme- 
try with rotation around the z-axis, but vv = and Vr-Q 



_ I dp <90 

r p dr dr 

1 dp 50 
- -TT" + ^ 

p oz oz 



(1) 
(2) 



where Vr, v. and are the three components of the velocity 
field, p the gas pressure and p the mass density, respectively. 
r - (x^ + y^y^^ is the distance from the symmetry axis and z 
is the distance from the midplane. Neglecting self-gravity, the 
gravitational potential is given by 



<l)(r, z) = - 



(3) 



where M* is the stellar mass and G the gravitational constant. 
We follow the ideaof "1-h1D" modelling (D'Alessio et al., 1998; 
Malbet et al., 2001 ; Duflemond et al., 2002) by assuming that the 
radial pressure gradient l/p{dp/dr) in Eq. (1) is small compared 
to centrifugal acceleration and gravity, in which case the radial 
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and the vertical components of the equation of motion decouple 
from each other. The radial component then simply results in 
circular Keplerian orbits 



xl/2 



(4) 



leaving the radial distribution of matter undetermined, as it is 
in fact mostly determined by the actual distribution of angular 
momentum in the disk. Consequently, the vertical component of 
the equation of motion (Eq. 2) can be solved independently for 
every vertical column in the disk 



Idp 
p dz 



zGM^ 



3/2 



(5) 



Equation (5) is integrated from the midplane upwards by sub- 
stituting the density for the pressure via p-c^p, and assuming 
that the isothermal sound speed cj is a known function of z. 
Numerically, we perform this integration by means of an ordi- 
nary differential equation (ODE) solver, using a simple point- 
wise linear interpolation of c^(z) between calculated grid points 
c\{rj, Zk)- Since for known cjiz) the solution p{z) has a free fac- 
tor, we put p{0) - 1 and scale the results later to achieve any 
desired colunm density at distance r 



E(r) = 2 



z)dz , 



(6) 



where the factor 2 is because of the lower half of the disk, which 
is assumed to be symmetric. In this paper, we assume a powerlaw 
distribution of the column density 



2(0 = So r- 



(7) 



in the main part of the disk, except for the "soft edges" (see 
Sect. 3.1), and determine 2o from the specified disk mass M^isk 



-2,1 



i;(r) r dr , 



(8) 



where Ri„ is the inner radius and Rom the outer radius of the disk. 
In summary, supposed that c^(r, z) is known, the disk structure 
is determined by the parameters Mdisk, A/*, Ri^, Rouu and e. 



(+ gravity). According to (Eq. 1) the force equihbrium in this re- 
laxed state is given by 



mRinY 



1 dp 50 

p dr dr 



(9) 



which provides an equation for the desired density structure p(r). 
Using p{r) = Cj p{r) and assuming = const, the result is 



In 



Pirp) 

"p(/?in) 



1 



(10) 



where vq is an arbitrary point inside R^. Generalizing Eq. (10) to 
column densities (with ci. measured in the midplane) we write 



2(ro) ~ 2(/?i„) exp 



1 



2r2 



(11) 



A similar expression can be found for the column density out- 
side of the outer boundary. The CO observations of Hughes 
et al. (2008) show that such treatments can improve model fits. 
However, we have chosen to apply our approach for soft edges 
only to the inner boundary in this paper. 

To summarize, if angular momentum is transported inside- 
out in the disk, the density structure may decrease more grad- 
ually or even increase further inward (Hartmann et al., 1998). 
However, it is hard to figure out any circumstances where the 
column density could decrease more rapidly at the inner rim as 
compared to Eq. (11). 

4. Continuum radiative transfer 

The chemistry and the heating & cooling balance of the gas in 
the disk (see Sects. 5 and 6) depend on the local continuous ra- 
diation field Jy{r, z) and the local dust temperature T^ir, z) which 
is a result thereof. These dependencies include 

1 . thermal accommodation between gas and dust, which is usu- 
ally the dominant heating/cooling process for the gas in the 
midplane (^ Tj), 

2. photo-ionization and photo-dissociation of molecules, as 
well as heating by absorption of UV photons, e. g. photo- 
electric heating (— > /uv), 

3. radiative pumping of atoms and molecules by continuum 
radiation which alters the non-LTE population and cooling 
rates, sometimes turning cooling into heating (— >/y), 

4. surface chemistry on grains, in particular the H2-formation, 
and ice formation and desorption (— > 7^, J\jw)- 



3.1. Soft Edges 

The application of a radial surface density powerlaw (Eq. 7) in 

the disk between Ri^ and 7?out is, although widely used, obviously 
quite artificial and even unphysical. Equation (1) demonstrates 
that an abrupt radial cutoff would produce an infinite force be- 
cause of the radial pressure gradient dp/dr, which would push 
gas inward at and outward at /?out, respectively, causing a 
smoothing of the radial density structure at the boundaries. 

Let us consider an abrupt cutoff in the beginning and study 
the motion of the gas as it is pushed inward due to the radial 
pressure gradient at the inner boundary. Since the specific angu- 
lar momentum L^iRm) =Rin v^iRm) = r v^(r) is conserved during 
this motion, the gas will spin up as it is pushed inward, until the 
increased centrifugal force balances the radial pressure gradient 



Previous chemical models have often treated these couplings 
by means of simplifying assumptions and approximate formula 
(e.g. Kamp & Bertoldi, 2000; Hollenbach et al., 1991; Nomura 
& Millar, 2005). For a rigorous solution, a full 2D continuum ra- 
diative transfer must be carried out, which provides 7d(r, z) and 
Jyir, z), including the UV part, at every location in the disk. 

ProDiMo solves the 2D dust continuum radiative trans- 
fer of irradiated disks by means of a simple, ray-based, long- 
characteristic, accelerated A-iteration method. From each grid 
point in the disk, a number of rays (typically about 100) are 
traced backward along the photon propagation direction, while 
solving the radiative transfer equation 

^-Sy-K (12) 
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assuming LTE and coherent isotropic scattering 
_ <''^Bv(7d) + CVv 

ly is the spectral intensity, Jv - j;; j U- dO. the mean intensity, 
5y the source function. By the Planck function, and i^^'^, /cj'^^ and 
/fj"' = K^'^ + Kf^ [cm"'] are the dust absorption, scattering and 
extinction coefficients, respectively. 

The dust grains of various sizes at a certain location in the 
disk are assumed to have a unique temperature 7^ in modified 
radiative equilibrium 



■ dust 



J Cjydv = J CByiTa)dv, 



(14) 



where the additional heating rate Fdust accounts for non-radiative 
heating (negative for coohng) processes hke thermal accommo- 
dation with gas particles and frictional heating. An accelerated 
A scheme is used to get converged results concerning Jv(r,z) 
and 7d(r,z). The details of this method will be described in the 
following sections. 

4.1. Geometry of rays 

Let ro = (xo, yo, zq) denote a point in the disk where the mean in- 
tensities Jv(ro) are to be calculated. The direction of a ray start- 
ing from ro is specified by a unit vector which points in the re- 
verse direction of the photon propagation 







' sin COS ' 








sin sin 


(15) 


. "3 . 




cos 6 





as specified in a local coordinate system where (0,0, 1) points 
toward the star. One ray is reserved for the solid angle occupied 
by the star as seen from point ro, subsequently called the 
"core ray". All other IxJ rays represent the remainder of the An 
solid angle by a 2D-mesh of angular grid points {0,- \ i = 1, /} 
and{0,|;= 1,...,7) 



Oi = e^ + (n-e^)\ — ] (16) 

<Pj = ^J^' (17) 

where 0* = arcsin(7?*/c/) is the half angular diameter of the star 
as seen from point ro,/?* the stellar radius, and (i = (Xy+>'g+ZQ)'''^ 
the radial distance, (p only ranges from to n, because the disk 
problem is symmetric /v(ro, 9, +<p)=Iv{rQ, 0, -<p). A power index 
p>\ (/? ~ 1 .5) assures that there are more rays pointing toward 
the hot inner regions than toward the cooler interstellar side. The 
integration over solid angle is carried out as 



!=1 ;=1 

Q* = 2;r(l - cose*) 



do., 



(pj) (cose,- -cos ^i+i) , 



(18) 

(19) 
(20) 



2(<^7+l 

The central direction of solid angle interval dD.ij is given by 0,- = 
(^i+i -I- 6i)/2 and = ((pj+i + 0j)/2, and these are the angles 
actually considered in Eq. ( 1 5). In order let the core ray with = 
point toward the star, we apply the following rotation matrix 



nx 

tly 



1 
, sina cos a 



(21) 



where a = y8 -i- 90" and tanyS = zo/(j(^ -i- j'^)'/^. 



4.2. Solutiori of the radiative transfer equatiori 

From every grid point ro along each ray in direction n = 
(tlx, riy, n.) we solve the radiative transfer equation (Eq. 12) back- 
ward to the photon propagation direction. The optical depth 
along the ray is given by 



Ty{s)= f C\ro + s'n)ds' 
Jo 

The formal solution of the transfer equation Eq. (12) is 

as)Syis)e-''^'^ ds 





(22) 



(23) 



where 7™ is the intensity incident from the end of the ray at Smax- 
We start a ray at * = with Ty = and ly = and choose a suit- 
able spatial step size As. For each step, the opacities and source 
functions (all wavelengths) at the start point and the end point 
of the step, ro-i-sn and ro + (s + As)n, are interpolated from the 
pre-calculated values on the grid points using a 2D-interpolation 
in cyhnder coordinates ((x^ + y^)^^'^, \z\). For the numerical in- 
tegration of Eqs. (22) and (23), we assume /c^"' = A + Bs and 
Sy-C exp(Ds), where the coefficients A, B, C, D are determined 
by the start and end point values. Simplifying the exponent by 
putting « A + BAs/2 yields 



+ iA + Bs')ds' 



Jo 



(24) 



/-A.S 

Iy(s+As) = Iy{s) + Ce-^-<'> (A + Bs') et^-^v^V _ (25) 

Jo 

The numerical integration is carried out with analytic expres- 
sions for these integrals. The procedure is repeated for two half 
steps of size Asf2. If the results differ too much, the step size As 
is reduced and the step is re-calculated. In case of small difl'er- 
ences, the step size is increased for the following step. 

At the end of each ray, the attenuated incident intensities 
y^mcg-TvCw) Qj-Q added according to Eq. (23), where for the core 
ray the stellar intensity P"'^ - I* is used, and for all other rays 
the interstellar intensity 7™ = 7'^"^ is applied. Non-core rays 
may temporarily leave the disk, but re-enter the disk after some 
large distance. These "passages" are treated with large, exactiy 
calculated As and zero opacities. 

For the 2D-interpolation, it turned out to be important to use 
a log-interpolation for the source function Sv{r,z) which can 
change by orders of magnitude, e. g. across a shadow, within 
one step. In case of linear interpolation, the numerical radiative 
transfer shows much more numerical difirision. 



4.3. Irradiation 

The radiation field in (and around) passive disks is completely 
determined by the stellar and interstellar irradiation, and the ge- 
ometry of the dust opacity structure. Therefore, setting the irra- 
diation as realistic as possible is of prime importance. 



Stellar irradiation For the incident stellar irradiation, a model 
spectrum from stellar atmosphere codes is used, e. g. a Phoenix- 
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T= 100000 



> 



1000 lO" 

X [nm] 

Fig. 2. Incident stellar intensity compiled from two sources: a Phoenix 
solar model spectrum with T^g = 5800 K, log g = 4.5, Z = 1 (black 
line) and the chromospheric flux of "young sun" HD 129333 (Dorren & 
Guinan, 1994, red line). Note that the ionizing and photo-dissociating 
flux is assumed to be restricted to the interval [91.2 nm, 205 nm] which 
includes Lycf at 121.6 nm. 

model'. Neglecting limb-darkening, the incident stellar intensi- 
ties are related to the surface flux at the stellar radius via 



(26) 



where F* is the spectral flux and Hy the Eddington flux. Young 
stars, which are active and possibly accreting, have excess UV as 
compared to model atmospheres, in particular cool stars. This is 
of central importance for ProDiMo , because it is just this radia- 
tion that ionizes and photo-dissociates the atoms and molecules 
in the disk. Therefore, we add extra UV-flux as e. g. reduced 
from observations or given by other recipes, see Fig. 2. 

Interstellar irradiation Assuming an isotropic intersteUar radia- 
tion field, all incident intensities for non-core rays are approx- 
imated by a highly diluted 20000 K-black-body field plus the 
2.7 K-cosmic background. 



rlSM 



ISM 



'■X 



1 .71 ■ Wdii Bv(20000 K) + By(2.7 K) 



(27) 



The applied dilution factor Wdii - 9.85357 x 10"'^ is calculated 
from the normalization - 1 according to Eq. (41), which is 
close to the value given by Draine & Bertoldi (1996). is a 
free parameter which describes the strength of the UV field with 
respect to standard interstellar conditions. 



4.4. Iteration and dust temperature determination 

In order to solve the condition of the dust radiative equilibrium 
(Eq. 14) and the scattering problem, a simple A-type iteration 
is applied. The source functions are pre-calculated on the grid 
points according to Eq. (13), with Jy - 7°''' and - 7'^''', and 
fixed during one iteration step. After having solved all rays from 
all points for all frequencies, the mean intensities are updated as 



1 

An 



/-I y-i 



/v(ro, 0, 0) Q* + ^ ^ /v(ro, 0,-, 0,) d^i^ 



'■=1 j=i 



, (28) 



and the dust temperatures are renewed according to Eq. (14). If 
the maximum relative change \ Jy - jf'^\l{Jy + 7smaii) (all points. 
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Fig. 3. Benchmark for the dust continuum radiative transfer part. 
Vertical cuts of the calculated dust temperature structure Ti(r,z) are 
shown at different radii of a disk with midplane optical depth Ti^^ = 10' 
at 1 yum. The ProDiMo results are shown with blue diamonds in com- 
parison to the MCFOST results (Pinte et al., 2006) with black lines. 



all frequencies, 7smaii = lO^-^'ergcm^^s^'Hz^'sr"') is larger than 
some threshold (~0.01), the source functions are re-calculated 
and the radiative transfer is solved again. In order to acceler- 
ate the convergence, we apply the procedure of Auer (1984). 
We benchmarked results of our radiative transfer method against 
results of other Monte-Carlo and ray-based methods in (Pinte 
et al., 2009). The convergence in optically thick disks is tricky, 
but we can manage test problems up to a midplane optical depths 
of about T- 10^ with this code (see Fig 3). 



4.5. Spectral bands and band-mean quantities 

The main purpose of the continuum radiative transfer in 
ProDiMo is to calculate certain frequency integrals, e. g. solv- 
ing the condition of radiative equilibrium for the dust grains 
(Eq. 12) or calculating the local strength of the UV radiation field 
X (Eq. 41). The incident stellar spectrum is strongly varying in 
frequency space, especially in the blue and UV (see Fig. 2) and 
the evaluation of these integrals, in principle, requires a large 
number of frequency grid points v^, which is computationally 
expensive. 

However, the incident radiation interacts with quite smooth 
and often completely flat dust opacities in the disk. Thus, it 
makes sense to "interchange" the order of radiative transfer 
and frequency integration, and to switch from a monochromatic 
treatment to a treatment with spectral bands. 

We consider a coarse grid of frequency points {vk\k = 
Q,...,K} (e.g. K = 12) which covers the whole SED, ranging 
from 100 nm to 1000 /im. Instead of By(T) we consider band 
means as 



see ftp://ftp.hs.uni-hamburg.de/pub/outgoing/phoenix/GAIA/ 



Bk(T) 



1 n 

— By{T)dV 

Jy,_, 



(29) 
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where Avj: = Vk-Vk-i-lna similar way, we treat the intensities, 
mean intensities and opacities 



4 



-f 

1 p 



Jvdv 



dv 



(30) 
(31) 
(32) 



Henceforth, we exchange the index y by the index k in all 
equations in this section, and retrieve the recipes for the band- 
averaged continuum radiative transfer. This is of course not an 
exact treatment, because it ignores all non-hnear couphngs, but 
an approximation that allows us to use fewer frequency grid 
points without loosing too much accuracy. 

4.6. Dust kind, abundance, size distribution, and opacities 

We assume a uniform dust abundance and size distribution 
throughout the disk. The dust particle density is given by 



rid 



I f(a)da, 

v/flniin 



(33) 



where a is the particle radius and f{a) is the dust size distribution 
function [cm which is assumed to be given by a powerlaw as 
/(a) = /consta • The moments of the size distribution are 



1 n^max 



f(d)da . 



(34) 



The constant in the powerlaw size distribution /const is deter- 
mined by the requirement that the dust mass density 



Pd = ndPgr—{a ) 



is given by a specified fraction of the gas mas density pd/p- pgr 
is the dust material mass density. 

The dust opacities are calculated from effective medium the- 
ory (Bruggeman, 1935) and Mie theory (Mmx from S. Wolf, ac- 
cording to Voshchinnikov, 2002). Any uniform volume mix of 
solid materials with known optical constants can be used. The 
dust opacities are calculated as 



na Qextia, X) fia) da , 



(36) 



where Qext(o, A) is the extinction efficiency. Similar formula ap- 
ply for absorption and scattering opacities, and k'^'^, where 
Qextia, A) is replaced by Q^hsia, ^) and gsca(a, /I), respectively. 

5. Chemistry 

The chemistry part of ProDiMo is written in a modular form 
that makes it possible to consider any selection of elements 
and chemical species. In the models presented in this paper, we 

consider chemical reactions involving A^ei = 9 elements among 
A^sp = 71 atomic, ionic, molecular and ice species as listed to 
Table 1. 

The rate coefiicients R are mostly taken from the Umist 2006 
data compilation (Woodall et al., 2007). Among the species 
listed in Table 1 we find 911 Umist chemical reactions, 21 of 
them have multiple Jg-fits. We add 39 further reactions which 



Table 1. Elements and chemical species 



9 elements 



71 species 



H, He, C, N, O, Mg, Si, S, Fe 



H, H+, H-, He, He+, C, C+, O, 0+, S, S+, Si, Si+, 
Mg, Mg+, Fe, Fe+, N, N+, H2, HJ, H*, H+, OH, 
0H+, H3O+, H2O, H2O+, CO, CO+, HCO, H2CO, 
HCO+, O2, OJ, CO2, C0+, CH, CH+, CH2, CHJ, 
CH3, CHJ, CH4, CH+, CHJ, SiO, SiO+, SiH, SiH+, 
SiH+, SiOH+, NH, NH+, NH2, NHJ, NH3, NH+, 
N2, HNJ, CN, CN+, HCN, HCN+, NO, N0+, C0#, 
H20#, C02#, CH4#, NH3# 



Ice species are denoted with "#". 

HJ designates vibrationally excited H2. 



are either not included in Umist or are treated in a more sophis- 
ticated way, as explained in Sects. 5.2 to 5.5. Among the alto- 
gether 950 reactions, there are 74 photo reactions, 177 neutral- 
neutral and 299 ion-neutral reactions, 209 charge-exchange re- 
actions, 46 cosmic ray and cosmic ray particle induced photo 
reactions, and 26 three-body reactions. The net formation rate of 
a chemical species i is calculated as 



jkt it 
V jkt jk 



(37) 



where Rjk^ie designates (two-body) gas phase reactions between 
two reactants j and k, forming two products i and {. Rf^^ji^ indi- 
cates a photo-reaction which depend on the local strength of the 
UV radiation field, and RjL,j( a cosmic ray induced reaction. 

5.1. Photo-reactions 



(35) Photon induced reaction rates can generally be written as 



(38) 



= 47: f a{v)^dv = 7 r (r{X)AuxdX 
J hv hj 

where cr(v) is the photo cross-section of the reaction, v the 
frequency, A the wavelength. It the Planck constant, and ua = 
^Ja the spectral photon energy density [erg cm "^J, respectively. 
The proper calculation of the photo-rates R^^ [s~'] according to 
Eq. (38) would require the calculation of a detailed (i. e. UV- 
line resolved) radiative transfer including molecular opacities to 
account for self-shielding eff'ects (to get Jy) as well as detailed 
knowledge about the wavelength-dependent cross section o-{A), 
which is not always available. 

In this paper, we will apply the Umist 2006 photo reaction 
rates in combination with molecular self- shielding factors from 
the literature instead. The application of detailed molecular UV 
cross sections in the calculated UV radiation field will be ad- 
dressed in a future paper. 

In the Umist database, photo-rates are given as 



,ph 



2x0 a ^Wi-yATn, 



(39) 



where xo is the unattenuated strength of the UV radiation field 
with respect to a standard interstellar radiation field, a the photo- 
rate in this standard ISM radiation field, and A™'^^ the extinction 
at visual wavelengths toward the UV light source. The photo- 
rates are derived for semi-infinite slab geometry, that is radiation 
is coming only from In; this explains the factor 2 in Eq. (39). 
In arbitrary radiation fields, for other than ID slab geometries. 
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r |AUJ r lAU] r [AUJ 



Fig. 4. Comparison of tiie UV radiation field strengths x (Eq. 41) between the simplified geometrical approach (l.h.s.), taking into account only 
vertical and radial dust extinction (see Eq. 47), and the results of the full 2D radiative transfer including scattering (middle) for the Mi;s]^ = 10"" Mq 
model. The two dashed contour lines show Ay = I and Ay = 10 as calculated from the minimum of the radial and vertical column densities. The 
r.h.s. figure shows the ratio ^/^'^''°, indicating that the regions mostly affected by scattering are situated roughly between Ay = 1 and Ay = 10. 



and for dust properties different from the ones used in UMIST, 
it is not obvious how to apply Eq. (39). In the following, we will 
therefore carefully explain our assumptions for the application 
of Eq. (39) to protoplanetary disks. 

Rollig et al. (2007) relate x - ^ a "unit Draine field" and 
we will follow this idea in ProDiMo. From the original work by 
Draine (1978), Draine & Bertoldi (1996) deduced 



4X 10" 



31.0164- 



49.913^., + 19.897 



4 



(40) 



for the standard ISM UV radiation field [erg cm ^] where 
/l/lOOnm. We apply an integral definition of^ as 



/-205r 



Aui dA 



1^205 r 
Jsi.lni 



Au^"''"'dA. 



(41) 



The wavelength interval boundaries have been chosen to en- 
sure coverage of the most important photo-ionization and photo- 
dissociation processes (van Dishoeck et al., 2006). Numerical 

integration yields Forame = i X'"',™ = 1-921 x 

1 0^^ cm"^ s" ' . Adopting the wavelength boundaries 91.2 nm and 
205 nm for the definition of our spectral band 1, we can directly 
calculate x from our banded radiative transfer method, including 
scattering, by 



47r , 

X ^ — Jl^Vl I -Foraine : 

nvi 



(42) 



where A^<h> is the hydrogen nuclei column density toward the 
UV light source and ^j"' the dust extinction coefficient per H- 
nucleus, averaged over spectral band 1. The ratio Ay/ruv de- 
pends on the frequency-dependence of the dust opacities as 



Av/t\j\ - 2.5 loge ■ K 



550iim''^l 



/4 



(45) 



However, in Eq. (39) we must not use Ay as calculated from our 
choice of dust properties ! Instead, we have to use the A™'^^ depth 
scale as used for the compilation of the Umist database. If we 
would use our Ay (from dust properties in the disk), the use of 
7 - containing Umist dust properties - would internally scale 
it to a Tuv that is wrong^. Thus, to obtain from our UV opti- 
cal depth the proper A™'^^, we need the ratio A™^^/tuv. Since 
the exact UMisT-ratio is not known, we calculate it according 
to Eq. (45) for standard "astronomical silicate" grains (Draine & 
Lee, 1984) with a size distribution /(a) oca"^-'' between 0.005 jim 
and 0.25 fim 



= 0.216 Tuv , 



(46) 



whereas for larger disk dust, a value around 1 is more typical. 

Another complicated problem is how to apply Eq. (39) in 
disk geometry. For this purpose we introduce a geometric mean 
intensity as it would be present, at least approximately, if only 
extinction but no scattering would occur 



(47) 



where we put vi - yvovi . The unshielded ISM photo-rate a is 
assumed to be given by 



a 



2hJ 



o-iA) Au^''"^"' dA 



(43) 



and the coefficient y in Eq. (39) can be identified as an effective, 
frequency-averaged opacity coefficient which contains implicit 
information about the frequency-range of the cross section cr(v), 
the Umist dust opacity and the shape of the ISM radiation field 
assumed. Neglecting gas extinction and assuming constant dust 
properties along the line of sight, the UV optical depth is 



Tuv 



<H> 



(44) 



/* and /p"^ are the incident band-mean stellar and interstellar 
intensities (see Sect. 4.3), and t^^ and rjjy are the radial (toward 
the star) and vertical (upwards) UV optical depth, respectively. 
Q* is the solid angle occupied by the star and Q}^^ -Att-Q.* the 
remainder of the full solid angle. Switching to corresponding x 
variables, we find 



n* 

ceo * — r^" 



ISM 



47r 



(48) 



' For species that can be ionized with visual light like H", it might be 
actually better to use our Ay scale. But most photo-reactions occur in 
the UV and Ay is just used as an auxiliary variable. 
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whereto is calculated with J\ - I* from Eq. (42), and;^^^ with 



Ji = ly ■ This decomposition into two slab geometries allows 
us to apply Eq. (39) and calculate the photo-rates as 

This approach to calculate the photo-rates according to Eq. (49) 
can be extended for molecular self- shielding factors (see 
Sect. 5 .2) and states a compromise between the usual two-stream 
approximation and a proper treatment of UV line-resolved radia- 
tive transfer according to Eq. (38). The factor xlx^''° corrects for 
the mayor shortcomings of the two-stream approximation, i. e. 
the effects of scattering and the assumptions about the geometry 
of the radiation field made in Eq. (47). Figure 4 shows that the 
enhancements xlx^'^° ^re very close to 1 in the upper, directly 
irradiated layers of the disk, but may be as large as 10^ at the in- 
ner rim and in the warm intermediate layer of the disk, and about 
10^ in the outer midplane due to scattering, or* is the unshielded 
photo-rate for 4;r-irradiation with star Ught, divided Xq- It can 
be calculated from cr{A) if known, or is assumed to be identical 
to 2a otherwise. 



5.2. Special UV photo reactions 

For the photo-ionization of neutral carbon, a = 3 x 10"'"s ' is 
taken from the Umist database, whereas a* is calculated accord- 
ing to the frequency-dependent stellar irradiation and the bound- 
free cross section of (Osterbrock, 1989). Molecular shielding by 
H2 and self-shielding is taken into account via the following fac- 
tors from Kamp & Bertoldi (2000) 



sc,c = Gxpi-cr^^Nc) 
,c.H. = exp(-0.9rr(^)"") 

■^c'' = *c,c sc,H2 Xo a exp(-Tuv) , 



(50) 
(51) 

(52) 



i. e. we refrain from an indirect formulation with Ay in cases 
we have the cross sections at hand. The approximation of H2 
shielding for the C ionization is strictly valid at low temperatures 
only (T < 300 K); for higher temperatures the factor 0.9 should 
be dropped. A^c and A^h, are the neutral carbon and molecular 
hydrogen column densities toward the UV light source, respec- 
tively, and cr^* = 1.1 X 10^'^cm^ the FUV-averaged cross section. 
The neutral carbon ionization rates due to radial and vertical UV 
irradiation are calculated separately according to Eq. (52), then 
multiplied by the respective solid angles and added together, and 
then corrected for scattering as in Eq. (49). 

For the photo-dissociation rate of molecular hydrogen, the 
same procedure applies with a H2 self-shielding factor taken 
from Draine & Bertoldi (1996, see their Eq.37). We assume 
Qr=4.2 X 10-"s-i (Draine & Bertoldi, 1996). 



■SH2,H2 - 



0.965 



0.035 



(1 -I- x/bn^y c exp(8.5 x lO"* c) 



(53) 



with X-NH2/5 X 10'"^cm~^, bu^ the H2 UV line broadening pa- 
rameter in [km/s] and c = (l +xy^^. The line broadening parame- 
ter b is defined as FWHM/(41n2)'^^; it contains the sum of ther- 
mal and turbulent velocities b = + Av^)'^^. Observations of 
line width in protoplanetary disks show that the turbulent veloc- 
ities are below 0.1 km/s (e.g. GuUloteau & Dutrey, 1998; Simon 
et al., 2000). 



The CO photo-dissociation rate is calculated from detailed 
band opacities in a similar fashion, taking into account the 
shielding by molecular hydrogen and the self- shielding. 



sco,H2 = exp(-5 ■ 10^) 



(54) 



with x = 9.555xlO-'^ (logipA^H^) -3.976. The CO photodisso- 
ciation rate for each band is interpolated from pre-tabulated rates 
using the CO column density, gas temperature and line broad- 
ening parameter as input (Kamp & Bertoldi, 2000; Bertoldi & 
Draine, 1996). 



5.3. H2 formation on grains 

The formation of H2 on grain surfaces H-i-H-i- grain Hi-n grain 
is taken into account according to (Cazaux & Tielens, 2002) 



= ^VniTs)nd^n{a^)aHeiTi) 



(55) 



with latest updates for the temperature-dependent efficiency 
e(Ti) from S. Cazaux (2008, priv.comm.). v* = {kTJ{2nmH)y'^ 
is the thermal relative velocity of the hydrogen atom, nj 4n{a^) 
is the total surface of the dust component per volume, and 
Oh ~ 0.223 is the sticking coefficient, which results in 
i v*(100K)4;r<fl2)(„^/„^jj^)a^ = 3 x lO^'^cm^s for standard 
ISM grain parameter pd/p - 0.01, - 0.005 /im, Omax = 
0.25 ;um, flpow = 3.5 and pgr - 2.5g/cm"^. The rate coefficient 
still needs to be multiplied by the neutral hydrogen particle 
density nn to get the H2 formation rate [cm~^s"']. 



5.4. Chemistry of excited H2 

13 reactions for vibrationally excited molecular hydrogen H* are 
taken into account as described in (Tielens & Hollenbach, 1985). 
The FUV pumping rate H2 + hv ^ H2* -H hv' is assumed to be 
10 times the H2 photo-dissociation rate. Two additional reactions 
are added for the coUisional excitation by H and H2 as inverse of 
the de-excitation reactions 

H2+H^H*-i-H: R = C^,(Tg) exp(-A£/A:Tg) (56) 

H2 -I- H2 ^ H* -I- H2 : R = C^liTg) expi-AE/kTg) , (57) 

where the energy of the pseudo vibrational level AE - 2.6 eV 
as well as the coUisional de-excitation rate coefficients C^^ and 
C^l [cm-^s-i] are given in (Tielens & Hollenbach, 1985). 



5.5. Ice formation and evaporation 

The formation of ice mantles on dust grains plays an impor- 
tant role for the chemistry in the dark and shielded midplane. 
At the moment, five ices are considered: C0#, C02#, H20#, 
CH4# and NH3# which are treated as additional species in the 
chemistry (Sects. 5 and 5.6). Apart from the adsorption and des- 
orption reactions of these species and the H2 formation on grains 
(Sect. 5.3) no other surface reactions are currently taken into ac- 
count. In particular, we do not form water on grain surfaces. 

Considering coUisional adsorption, and thermal, cosmic-ray 
and photo-desorption, the total formation rate of ice species i is 



dn 



'* „ Dads „desorb / ndes.th , r>des,ph „des,cr\ 

— = HiK^ - n.# \K. + K, + K 



(58) 



where n,# is the density of ice units / and nf^""'^^ the fraction of n,# 
located in the uppermost active surface layers of the ice mantle. 
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5.5.1. Adsorption 

A gas species will adsorb on grain surfaces upon collision. The 
adsorption rate [s"'] is the product of the sticking coefficient a, 
the total grain surface area per volume 4n{a^)nd and the thermal 
velocity v* = {kT^I{2nmd)^'^ 



■ An{a^)ndavf 



(59) 



where m/ is the mass of gas species ;. We assume unit sticking 
coefficient {a -I) for all species heavier than Helium (Burke & 
HoUenbach, 1983). 



Table 2. Adsorption energies and photo-desorption yields. 







T^lir»fr»-Hf»«nt*nti ATI viplH 




[K] 


y, [per UV photon] 


co# 


960" 


2.7 X 10-3 c 


C02# 


2000" 


1.0 X 10-3 " 




4800* 


1.3 X 10-3 


CH4# 


1100" 


1.0 X 10-3 « 


NH3# 


880" 


1.0 X 10-3 « 



a: Aikawa et al. (1996); b: Sandford & AUamandola (1990); c: 
Oberg et al. (2007); d: Westley et al. (1995); e: assumed values. 



5.5.2. Desorption 

A chemical species with internal energy greater ffian the energy 
that binds it to a grain surface will desorb. Desorption mecha- 
nisms depend on the source of the internal energy. 

I. Thermal desorption: An ice species ; at the surface of a grain 
at temperature Tj has probability to desorb 

^ads N 
/ 



R 



ides.th 



vf'^exp 



(60) 



where v^**^ = (2 kEf^ l{Ti^mi)yi^ is the vibrational frequency 
of the species in the surface potential well of ice species i, Wsurf = 
1.5 X 10^^ cm"^ is the surface density of adsorption sites and Fiinj) - 
ZiJ"^^ is the adsorption binding energy. The adopted values are 
provided in Table 2. Following Aikawa et al. (1996), ffie number 
density of ice units at the active surface n^*"* is given by 



where fcR is the cosmic ray ionization rate of H2, filOK) = 
3 . 1 6 X 1 0- the 'duty-cycle' of ffie grain at 70 K and ''*(70 K) 
ffie ffiermal desorption rate for species i at temperature = 70 K. 
The adopted value for /(lOK) is strictly valid only for 0.1 fxm 
grains in dense molecular clouds. 

5.6. Kinetic cliemical equilibrium 

Assuming kinetic chemical equilibrium in the gas phase, and be- 
tween gas and ice species, we have ^ = in Eq. (37) and obtain 
i = l ... Nsp non-linear equations for the unknown particle densi- 
ties Hj (7 = 1 ...A^sp) 



(64) 



It is noteworffiy that ffie electron density is not among ffie un- 
knowns, but is replaced by ffie constraint of charge conservation 



desorb 



ice ^ 
' '^tot - '^act 



(61) 



(65) 



where Wact - '\-n{a )nd Wsurf A^Lay is the number of active surface 
places in ffie ice mantle per volume and A^Lay is the number of 
surface layers to be considered as "active". We assume A^Lay = 2 
in accordance wiffi (Aikawa et al., 1996). nf^^ - liyn;# is ffie 
number density of ices. 

2. Photo-desorption: Absorption of a UV photon by a surface 
species can increase ffie species internal energy enough to induce 
desorption. The photo-desorption rate of species i is given by 



ndes.ph / 2\ ''d v 

=n{a} Yi xFuraine 

lact 



(62) 



where zj is the charge of species j in units of the elementary 
charge. The expUcit dependency of ric on ffie particle densities 
Hj causes additional entries in ffie chemical Jacobian dFi/dnj = 
dFj/drij + dFi/drie ■ dn^ldrij. 



5.7. Element conservation 

The system of Eqs. (37) is degenerate because every individ- 
ual chemical reaction obeys several element conservation con- 
straints, and ffierefore, certain Unear combinations of Fs can 



where 7, is the photo-desorption yield (see Table 2), x^Tytsme 
is a flux-like measure of the local UV energy density 
[photons/cm^/s] computed from continuum radiative transfer 
(Eqs. 41, 42). Photo-desorption can enhance gas-phase water 
abundances by orders of magnitude in outer region of disks n^H> ^k—/ i; Vi,k = 
(Willacy & Langer, 2000; Dominik et al., 2005; Oberg et al., i 
2008). 

3. Cosmic-ray induced desorption: Cosmic-rays hitting a grain 
can locally heat ffie surface and trigger desorption. Cosmic-rays 
can penetrate deep into obscured regions, maintaining a mini- 
mum amount of species in ffie gas-phase. Cosmic-ray fluxes in 
disks may be higher ffian in molecular clouds because of ffie 
stellar energetic particles in addition to the galactic component. 
X-ray photons can also penetrate deep inside ffie disk and locally 
heat a dust grain but X-ray induced desorption is not included in 
the code yet. We adopt for the cosmic-ray desorption the formal- 
ism of Hasegawa & Herbst (1993). 



be found which cancel out, making the equation system under- 
determined. Only if ffie element conservation is implemented in 
addition, ffie system (Eqs. 37) becomes well-defined. 

Considering the total hydrogen nuclei density n<H> as given, 
ffie conservation of element k can be written as 



(66) 



resulting in A'ei auxiliary conditions, ek is the elemental abun- 
dance of element k normalized to hydrogen and v, ,t are the stoi- 
chiometric coefficient of species i with respect to element k. 

Alternatively, the gas pressure p may be considered as ffie 
given quantity and ffie relative element conservation can be ex- 
pressed by 



hint We + ^ njm,- j - ^ v,-,t mk = . 



(67) 



R 



ides.cr 



fC70K)R'^^'''\70K) 



5 X lO-i^s-i 



(63) 



where q = ekmk/CLk' ^k'tn^') is the relative mass fraction of ele- 
ment k, nij the mass of a gas particle of kind / and mi^ the mass of 
element k. Since surmning up all Eqs. (67) for A: = 1 ... A^ei results 
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Table 3. Assumed element abundances in (gas + ice) 



element 


12 + log 6 


mass 11 action e 


H 


12.00 


7.66x10' 


He 


10.88 


2.28x10^' 


C 


8.11 


1.19x10"^ 


N 


7.33 


2.28x10-"* 


O 


8.46 


3.52x10-3 


Mg 


6.62 


7.76x10-5 


Si 


6.90 


1.70x10-* 


S 


6.28 


4.64x10-5 


Fe 


6.63 


1.82x10-'' 



This choice of element abundances implies p = 1.315 amu • n^ny- 



in p =p, one of these equations is redundant and can be replaced 
by the constraint of given pressure 



p- +z,)^rg = 0. 



(68) 



where k the Boltzmann constant. The element conservation 
is implemented by replacing A^ei selected components of Fj 
in Eq. (64) by these auxiliary conditions, either according to 
Eq. (66) or according to Eqs. (67) and (68), after suitable nor- 
malization. For this purpose, we choose for every element k the 
index j that belongs to the most abundant species containing this 
element. Eq. (68) overwrites the entry for the most abundant H- 
containing species. 

The global iteration, which solves the hydrostatic disk struc- 
ture consistently with the chemistry and heating & cooling bal- 
ance (see Fig. 1), is found to converge only if the chemistry 
is solved for constant pressure p. Since the vertical hydrostatic 
condition (Eq. 2) is a pressure constraint, it is essential to en- 
sure that the chemistry solver, coupled to the Tg -determination 
via heating & cooling balance, is not allowed to change p as it 
would be the case if n(H) was fixed. At given pressure p, Tg may 
be found to increase during the course of the iteration, but only 
if simultaneously n<H> drops, thereby conserving the /^-structure 
within one global iteration step. 



5.8. Numerical solution of chemistry 

The non-linear equation system (64), expressing the kinetic 
chemical equilibrium including element conservation, is usu- 
ally solved by means of a self-developed, globally convergent 
Newton-Raphson method. A quick and reliable numerical solu- 
tion of the Eqs. (64) is crucial for the computational time con- 
sumption, stability, and global convergence of our model. Our 
numerical experience shows that a careful storage of converged 
results (particle densities) is the key to increase stability and per- 
formance. These particle densities are used as initial guesses for 
the next time the Newton-Raphson method is invoked, either in 
form of a downward-outward sweep through the grid (first it- 
eration), or from the last results of the same point (following 
iterations). 

In cases, where the solution by the Newton-Raphson method 
fails, we fall back to the time-dependent case and solve Eqs. (37) 
by means of the ODE solver Limex (Deuflhard & Nowak, 1987) 
for 10^ yrs, which is much slower but in practice gives the same 
results as the Newton-Raphson method. 




Fig. 5. Different pumping and escape probabilities according to the pre- 
dominantly radial irradiation and the predominantly vertical escape. 



6. Gas thermal balance 

The net gain of thermal kinetic energy is written as 



dt 



(69) 



where Fjt and are the various heating and cooling rates 
[ergcm-^s-'] which are detailed in the forthcoming sub- 
sections. Restricting ourselves to the case of thermal energy bal- 
ance, we assume dejdt -Q in the following and Eq. (69) states 
an implicit equation for the unknown kinetic gas temperature 7g. 
Since the heating and cooling rates depend not only on 7g, but 
also on the particle densities «sp, which themselves depend on 
Tg, an iterative process is required during which T„ is varied and 
the the chemistry is re-solved until Ta satisfies Eq. (69). 



6.1. Non-LTE treatment of atoms, ions and molecules 

The most basic interaction between matter and radiation is the 
absorption and emission of line photons by a gas particle, which 
can be an atom, ion or molecule. We consider a A^-level system 
with bound-bound transitions only and calculate the level popu- 



lations nj [cm ] by means of the statistical equations 



(70) 



which are solved together with the equation for the conservation 
of the total particle density of the considered species 2, rii-n 
The rate coefficients are given by (Mihalas, 1978): 



sp- 



Rul - Aul + BulJul + Cul 
Rlu = BluJul + Clu , 



(71) 



where u and / label an upper and lower level, respectively. A„/, 
Bui, B/u, Cul and C;„ are the Einstein coefficients for spontaneous 
emission, absorption, stimulated emission and the rate coeffi- 
cients for collisional (de-)excitation, respectively. Additionally 
we have the Einstein relations B„//A„/ = c^/(2/iv^p, Biu/Bui = 
g„/gi and the detailed balance relation C;„/C„z = gu/gi ■ 
exp{-AE,ii/kTg), where v,,/, gu, gi and AEui are the line center 
frequency, the statistical weights of the upper and lower level 
and the energy difference, respectively. The line integrated mean 
intensity is given by 

7„/ My, n) I An) dv dn (72) 

where <p„i{v, n) is the line profile function in direction n. 

6.1 .1 . Escape probability treatment 

The spectral intensity ly in Eq. (72) is affected by line absorption 
and emission. Assuming that the line source function (Eq. 74) 
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varies slowly in a local environment where the line optical 
depths (Eq. 75) build up rapidly, we can approximate for a static, 
plane-parallel medium 

. CV)exp(^)..^l-exp(^::^^))(73) 

where /JJ^Xju) is the continuous background intensity which 
propagates backward along the ray in direction fi. The direction 
//=1 points "outward" (/i = cos0). The line source function and 
the perpendicular line optical depth are given by (Mihalas, 1978) 



id 



2hv 



ul 



- 1 



SttvIAv 

It! 



D Jz ^ 81 ' 



(74) 
(75) 



where Av^ is the (turbulent + thermal) velocity Doppler width 
of the line, assumed to be constant along the line of sight in 
Eq. (75). Equations (72) and (73) can be combined to find 



ppump jco] 



esc 
Pul 



(76) 
(77) 




1000.0 



Fig. 6. Continuum mean intensities as input for non-LTE modelling. 
The calculated band-mean mean intensities are shown for one particular 
point (r,z) in a model (12 black dots) and a cubic spline interpolation 
through these points (black line). The vertical lines indicate the interval 
boundaries of the 12 spectral bands. The red line shows the band-mean 
incident stellar intensities vl* and the blue line shows the incident inter- 
stellar intensities v/'^*^. The radiation field has two major components, 
the dust attenuated UV - near IR part, originating mainly from the star, 
and the self-generated mid - far IR part, originating from thermal dust 
emission in the disk. 



where the direction-dependent and the mean escape probabilities 
are found to be 



pesc 



4>(x) exp ( - 

oo 



) dx 



(78) 
(79) 



with dimensionless line profile function (p(x) = exp(-x^)/y/7T, 
X = {v - Vui)/Av£, and frequency width Avq - v„/Av£)/c. Using 
Eq. (77), it is straightforward to show that the unknown line 
source function S^^j can be eliminated, and the leading term 
riuAui cancels out, when considering the net rate n„Au; -i- (riuBui - 
niBi,)J,, = n,A,iP^f' + (n.B,, - n,B,„)pP™V™"'. Thus, we can 
solve the statistical rate Eqs. (70) with modified rate coefficients 



5 _ D pPumP rcont , ^ 
Kill - "lu^^iii Jv„i + ^hi , 



(80) 



which is known as escape probability formalism (Avrett & 
Hummer, 1965; Mihalas, 1978). P is the mean probability for 
line photons emitted from the current position to escape the lo- 
cal environment and p the mean probability for continuum 
photons to arrive at the current position. 7^;°"' is the continuous 
mean intensity at line center frequency v„/, as would be present 
if no line transfer effects took place. In semi-infinite slab sym- 
metry, all directions yU < have infinite line optical depth and can 
be discarded from the calculation of the escape probabilities 



P:r(r7) = X 0W exp{-^ii^)di,dx 

^ J-oo Jo P 

= - <t>ix)E2{rlY<p(x))dx 
This function is numerically fitted as 



(81) 
(82) 



0.5 ,T<0 
0.5 + (0.1995 lnT-0.2484)T- 0.04594t2 ,0<t<0.6 
1 -exp(-1.422T) 



3.324t + 0.2852t2 
0.1999/ 



-(ln(0.4799T)) 



-0.4195 



,0.6<T<9 



,9<T 



Considering the pumping probability as defined by Eq. (77), it 
is noteworthy that 7'^™'' ~ 7',"'^ is only valid in an almost iso- 
tropic background radiation field. In disk symmetry, much of the 
pumping is due to direct star light (see Fig. 4) which has a very 
pointed character. In the optically thick midplane, the contin- 
uous radiation field is almost isotropic, but here the pumping is 
pointless, because the radiation is thermalized and the collisional 
processes dominate. Considering near to far IR wavelengths at a 
certain height above the midplane, the irradiation from under- 
neath plays a role, but these directions are just the opposite of 
what is considered in Eq. (82), and so using 7',^™'' 
be strongly misleading. Thus, we approximate 



would 

III 



X+co 
<P(x)cXp{-T^f<P(x))dx 
OO 



(83) 



with t"'* now being the radially inward line optical depth. This 
function is numerically fitted as 

1 ,r<0 
1 -0.3989t-^0.09189t2-0.01497t3 ,0<t<0.9 
1 - exp(-0.6437T) 



0.6295t-h0.07008t2 
(^(ln(0.3367r)p°^ 



,0.9<T<9 
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6.1.2. Background radiation field 

The continuum background mean intensities have an im- 
portant impact on the gas energy balance. For example, in strong 
continuum radiation fields, the reverse process to line cooling, 
namely line absorption followed by collisional de-excitation, 
dominates. is identified to be just given by the mean in- 
tensities calculated from the dust continuum radiative transfer 
(see Sect. 4). In order to obtain the required monochromatic 
mean continuum intensities at the line center positions, we ap- 
ply a cubic spline interpolation to the calculated local continuum 
Jy°'^\r, z) in frequency space as depicted in Fig. 6. 



6.1 .3. Solving the statistical equations 

Equations (70, 75, 80) form a system of coupled equations 
for the unknown population numbers n,. Since the fine optical 
depths (Eq.75) depend partly on the local n„ these equations 
must be solved iteratively. We apply a fully implicit integration 
scheme for the numerical solution of (Eq. 75) where the line op- 
tical depth increment along the last downward integration step, 
between the previous and the current grid point, is given by the 
local populations, i. e. 



Table 4. non-LTE model atoms, ions and molecules 



8;n/3;AvB 



«<H> 



(^<H>(z,)-^(H>(z,-+i)). (84) 



where A^<h>(z!) is the vertical hydrogen column density at grid 
point Zi, and n„, «/ and «<h> refer to the current grid point i. A 
simple A-type iteration scheme is found to converge within typ- 
ically to 20 iterations. The outward radial line optical depth 
increments At™'^ between rj-\ and rj are calculated in a similar 
fully implicit fashion. 

6.1 .4. Calculation of the heating/cooling rate 

Once the statistical equations (Eqs. 70) have been solved, the ra- 
diative heating and cooling rates can be determined. There are 
two valid approaches. For the net cooling rate, one can either 
calculate the net creation rate of photon energy (radiative ap- 
proach), or one can calculate the total destruction rate of thermal 
energy (collisional approach). 



A, 



rrad = Yj"'^^"'^"' 

u>l 

U>1 
U>1 

niCiuAEui 



pump 
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(85) 



(86) 



(87) 



(88) 



U>1 



Both approaches must yield the same net result Frad-Arad = Fcou- 
Acou, which can be used to check the quahty of the numerical 
solution. In practise, one pair of these heating/cooling rates is 
often huge in comparison to the other pair, e. g. {Frad, Arad) » 
{Fcoi, Acoi} in a thin gas with radiatively controlled populations, 
and {Frad, Arad) ^ {Fcol, Acol) in a dense gas with close to LTE 
populations. Thus, it is numerically favorable to choose 



F 
A 



Fcol 
Trad 

Acol 

Arad 



> Trad > Fcol 
, Trad ^ Tcol 

, Frad > Fcol 
, Frad ^ Tcol 



(89) 
(90) 



species 


#levels #Iines 


coll. partners 


reference 


01 


3 


3 


p-H2,o-H2,H,H+,e- 


A-database 


CI 


2 


1 


p-H2,o-H2,H,H+,He,e- 


A-database 


CII 


3 


3 


H2,H,e- 


A-database 


Mg II 


8 


12 


e" 


Chianti 


Fell 


80 


477 


e" 


Chianti 


Sill 


15 


35 


e" 


Chianti 


sn 


5 


9 


e" 


CfflANn 


CO rot. & ro-vib. 


110 


243 


p-H2,o-H2,H,He,e- 


see text 


0-H2 ro-vib. 


80 


803 


p-H2,o-H2,H,He 


see text 


P-H2 ro-vib. 


80 


736 


p-H2,o-H2,H,He 


see text 


0-H2O rot. 


45 


158 


H2 


A-database 


P-H2O rot. 


45 


157 


H2 


A-database 



6.1.5. Atomic and molecular data 

The atomic and molecular data for Oi, Ci, Cn and H2O (en- 
ergy levels Ei, statistical weights gi, Einstein coefficients Aui, 
and collision rates C,,/ are taken from Leiden's LAMBDA-database 
(Schoier et al., 2005), see Table 4. In addition to these low- 
temperature coolants, we have included several ions as high- 
temperature coolants from the CfflANTi-database (Dere et al., 
1997): Mgii, Fen, Sin and Sn, taking into account all energy 
levels up to about 60 000 cm"' . This database has collisional data 
for free electrons only, but since we consider only ions of abun- 
dant elements here, the electron concentration is always rather 
high wherever these ions are abundant. Since the electron col- 
lisional rates are typically 10* times larger than those of heavy 
particles, the thereby introduced error seems acceptable. 

For CO, we have merged level and radiative data (£,, gi and 
Aui) of the rotational and ro-vibrational states (v = 0, 1,2, 3,4) 
from the Hitran database (Rothman et al., 2005) with collisional 
data among the rotational levels from the Lambda database. For 
the vibrational collisions we use the Ci^o data for H and H2 de- 
exciting colUsions from (Neufeld & Hollenbach, 1994) and for 
He collisions from (MilUkan & White, 1964). The de-exciting 
rate coefficients for other than 1^0 vibrational transitions are 
estimated according to the formula provided by (Elitzur, 1983) 



Cv'-»v = (v - v') Ci_,o exp 



(V 



1)1.5 eiT^ 



1 + \.5 6/Ta 



(91) 



where 6 = fko/k is the difference between the first vibrationally 
excited and the ground state energy. For detailed ro-vibrational 
modelling, these total vibrational collisional rates still need to 
be spread over the rotational sub-states. For simplicity, however, 
we assume Cvv^v/ ~ O-tv for every rotational state J. 

For H2, the level and radiative data (quadrupole transitions) 
are taken from Wolniewicz et al. (1998). We include calculated 
collisional excitations by H (Wrathmall et al., 2007), ortho- and 
para-H2, and HeUum (Le Bourlot et al., 1999). The H2 and H2O 
ortho to para abundance ratios are assumed to be at thermal equi- 
librium according to the gas temperature. 

6.2. Specific tieating processes 

Below, we list further heating processes that are not covered 
by Sect. 6.1. Photoelectric heating, cosmic ray ionization, car- 
bon photo-ionization and H2 photo-dissociation are stiU radia- 
tive processes, while other heating mechanisms are of chemical 
nature, such as H2 formation heating, or of dynamical nature, 
such as viscous heating. 
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6.2.1. Photo-electric heating 

UV photons impinging on dust grains can eject electrons with 
super-thermal velocities which then thermalize through colli- 
sions with the gas. The efficiency of this process decreases 
strongly with grain charge (positively charged grains are less ef- 
ficient heaters). The grain charge is set by the balance of incom- 
ing UV flux that ejects electrons and collisional recombination. 
The colUsion rate for recombination scales with electron density, 
thermal velocity and the ratio between potential and thermal en- 
ergy (O = ell/kTg, with U being the grain potential). Thus the 
grain charge can be parameterized by a 'so-called' grain charge 
parameter (Bakes & Tielens, 1994) 



(92) 



The probability of electron ejection after photon absorption 
(yield), is generally taken from experimental data on bulk mate- 
rial with large flat surfaces, and then applied to (smaller) astro- 
physical dust grains to compute the photoelectric heating rates. 
The heating process is thought to be less effective for micron- 
sized grains compared to small ISM dust grains. The reason 
is that a photo-electron can more easily be trapped within the 
matrix of a large grain, thus lowering the photoelectric yield. 
Experimental data on realistic astrophysical dust grains is sparse 
and only recently (Abbas et al., 2006) carried out experiments 
with sub-micron to micron sized individual dust grains. They 
measure yields that are larger than those of bulk flat surfaces 
and they find an increasing yield with increasing grain size. 
However, the underlying physics of these experiments are not 
yet well understood. 

Kamp & Bertoldi (2000) provide a formula to approxi- 
mate the photoelectric heating rate for large graphite and sili- 
cate grains using the photoelectric yields of bulk material from 
Feuerbacher & Fitton (1972). For silicate grains, the photoelec- 
tric heating rate FpE and the efficiency e are 



FPE = 2.5 X 10-V^%<H> e;r 



e = 



y = ' 



0.06 



1 -1-1.8 X 10-3x0-91 

0.7 ,x<10-'' 
0.36 , 10-4 <x<l 
0.15 ,x>l 



1-1- O.Olx 



(93) 
(94) 

(95) 

,-3 



valid for electron particle densities 10" cm" < We < 10 cm 
gas temperatures lOK < Tg < 10000 K, and strength of FUV 
radiation field 10"^ <X< 10^. Here, is the grain absorption 
cross section per H-nucleus (£r^''*n<H> =i^^, see Eqs. 32 and 36). 



6.2.2. PAH heating 

Very small dust grains such as polycyclic aromatic hydrocarbons 
(PAHs) are an extremely efficient heating source for the gas. 
The photoelectric heating rate can be written separately from the 
rest of the grain size distribution using only the first term of the 
(Bakes & Tielens, 1994) efficiency formulation 



FpE = /pah lO-^n(H) ex 
where the efficiency e is given by 
0.0487 



l+A-x 10-3x0-73 



(97) 



In the ISM, the abundance of PAHs is /pah = 1 ■ For disks, this 
value can be scaled according to the observed strength of the 
PAH bands. 

6.2.3. Carbon photo-ionization 

Ionization of carbon releases electrons with energies around 
1 eV (Black, 1987). Subsequent collisions heat the gas as 

Fc = 1.602 X IQ-^^Rfnc (98) 

where the photo-ionization rate R^^ is given by Eq.(52). 

6.2.4. Hi photo-dissociation heating 

Photo-dissociation of molecular hydrogen occurs via UV line 
absorption into an electronically excited state followed by spon- 
taneous decay into an unbound state of the two hydrogen atoms. 
The kinetic energy of such H-atoms is typically 0.4 eV (Stephens 
& Dalgarno, 1973), leading to an approximate heating rate of 

(99) 



r"? = 6.4 X 10 

ph 

Here 



-13 ^H. 

ph 



R^^ is the H2 photo-dissociation rate given in Sect. 5.2 
including dust and H2 self- shielding. 



6.2.5. cosmic ray heating 

Cosmic rays have a typical attenuation depth of 96gcm--^ 
and thus reach much deeper than stellar FUV photons (~ 

10~3gcm~^, see Bergin et al., 2007, for an overview). They 
ionize atomic and molecular hydrogen and this inputs approx- 
imately 3.5 and 8 eV into the gas for H and H2, respectively 
(Jonkheid et al., 2004). The heating rate can then be written as 



FcR = -TcR (5.5 X lO-'^ne + 2.5 x 10-"«h,) (100) 
where ^cr is the primary cosmic ray ionization rate. 

6.2.6. H2 formation heating 

The formation of H2 on dust surfaces releases the binding en- 
ergy of 4.48 eV. Due to the lack of laboratory data, we follow 
the approach by Black & Dalgamo (1976) and assume that this 
energy is equally distributed over translation, vibration and rota- 
tion. Hence, about 1 .5 e V per reaction is liberated as heat 



->H2 

form 



2.39 X 10^12 ^H2«H 



(101) 



where the H2 formation rate R\i^ is given in Sect. 5.3. 



6.2.7. Heating by collisional de-excitation of H* 

The fluorescent excitation of H2 by UV photons H2 -l- hy — > 
H** H* + hv' produces vibrationally excited molecular hy- 
drogen H* (Tielens & Hollenbach, 1985), and the vibrational 
excitation energy can be converted into thermal energy by colh- 
sions. The heating rate is 



'■ coll 



= ^ER 



coll 



"H* - "Hz exp 



( 



(96) <-.Hj = «H cS(rg) + nH, c«/(rg) 



(102) 
(103) 



where the excitation energy of the pseudo vibration level A£ 
and the collisional de-excitation rates C„; are given in (Tielens 
& Hollenbach, 1985), see also Sect. 5.4. The second term in 
Eq. (102) corrects for collisional excitation. 
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6.2.8. Viscous heating 

Due to high optical thickness, radiative heating cannot penetrate 
etiiciently to the midplane. These dense layers can instead also 
be heated by local viscous dissipation (Frank et al., 1992) 



Table 5. Parameter of the model depicted in Figs. 7 to 14. 



2 

rvis = ^P'^'kin^^kep 



(104) 



In the absence of a well-understood mechanism, angular mo- 
mentum transport is conceptualized using the kinematic - or 
turbulent - viscosity van often parameterized as an a- viscosity 
(Shakura& Syunyaev, 1973) 



Vkin = a cr H„ 



(105) 



where a is a dimensionless scaling factor, - pip is the isother- 
mal sound speed, Hg - cjl^^ef is the gas scale height, and 
f^kep = v^lr is the Keplerian angular velocity (see Eq.4). For 
r^2 AU, the viscous heating is known to be capable of dominat- 
ing the energy balance in the midplane (D' Alessio et al., 1998)^. 



quantity 


symbol 


value 


stellar mass 




IMo 


effective temperature 


Tts 


5770 K 


stellar luminosity 




lio 


disk mass 




0.01 Mo 


inner disk radius 


^in 


0.5 AU<" 


outer disk radius 


^out 


500 AU 


radial column density power index 


e 


1.5 


dust-to-gas mass ratio 


Pd/P 


0.01 


minimum dust particle radius 




0.1 pm 


maximum dust particle radius 


^max 


lOyum 


dust size distribution power index 




2.5 


dust material mass density 


Pgr 


2.5 gcm"^ 


strength of incident ISM UV 


^SM 


1 


cosmic ray ionization rate of H2 


fcR 


5 X 10-" S-' 


abundance of PAHs relative to ISM 


/pah 


0.12 


a viscosity parameter 


a 






(1) : soft inner edge applied, see 

(2) : dust optical constants from 



Sect. 3.1 

Draine & Lee (1984) 



6.3. Specific cooiing processes 

Most cooling processes are radiative in nature and covered in 
Sect. 6.1. However, two prominent high temperature cooling 
processes are treated in a simpler approximative fashion: Lyman- 
a and Oi-630nm cooling. 



6.3.1. Ly-a cooling 

Cooling through the Lyman-o; line becomes efficient at temper- 
atures of a few 1000 K (Sternberg & Dalgamo, 1989). Given the 
densities of atomic hydrogen riu and electrons n^, the cooling 
rate can be written as 



ALy-a = 7.3 X 10-iVneexp(-118400/rg) . 



(106) 



6.3.2. Oi-630nm cooling 

Line emission from the meta-stable level of neutral oxy- 
gen efficiently cools the gas at temperatures in excess of a few 
1000 K. With «o denoting the neutral oxygen particle density 
the cooling rate is (Sternberg & Dalgarno, 1989) 

Aoi630mn = 1-8 X 10^24^^ n, cxp (-22 800/rg) . (107) 

6.4. Misceiianeous tieating/cooiing processes 

We list below two additional processes that can cause either heat- 
ing or cooling of the gas. 



6.4.1. Thermal accommodation 

Following (Burke & HoUenbach, 1983), the energy exchange 
rate by inelastic collisions between grains and gas particles is 



Fg-g = 4 X 10"i2 n{a^) na n<H> orx V^g (^d - T^) 



(108) 



^ Without further adjustments, the viscous heating rate according to 
Eq. (104) scales as Tvis oc p at given radius r. Since all known cooling 
rates scale as A oc in the low density limit, there is always a critical 
height z above which the viscous heating would dominate the energy 
balance and lead to ever increasing Tg (well above 20000 K) with in- 
creasing height z- We consider this behavior as an artefact of the concept 
of viscous heating and/or a-viscosity. 



For gas temperatures Tg higher than dust temperatures Id. this 
rate turns into a cooUng rate Ag g. The thermal accoimnodation 
coefficient aj is set to the typical value for silicate and graphite 
dust of 0.3 (Burke & HoUenbach, 1983). 

6.4.2. free-free heating/cooling 

Free-free transitions directly convert photon energy into thermal 
energy (ff-heating) or vice versa (ff-cooling) during electron en- 
counters. The heating rate Fff and cooling rate Aff are given by 



Fff = 4;r 



: J K^Jydv (109) 
Aff = 4;r J ^.B^T^) dv (110) 

= nlcr^+ Me «H of' + Me "Hj 0\^'^ + «e "He of ^"'^ , (1 1 1) 

where is the free-free gas opacity [cm^^]. The free-free cross- 
sections (tJJ [cm^] for bremsstrahlung of singly ionized gases are 
taken from (Hummer, 1988), for H-fFfrom (Stilley & Callaway, 
1970), for H-fF from (Somerville, 1964), and for He-ff from 
(John, 1994). 

7. Sound Speeds 

After the chemistry (see Sect. 5) and the thermal gas energy bal- 
ance (see Sect. 6) have been solved throughout ffie disk volume, 
all particle densities and the kinetic temperature of the gas 
Tg are known, and ProDiMo can update the isothermal sound 
speeds on the numerical grid c^irj, Zk) as preparation for the next 
iteration of the hydrostatic disk structure (see Sect. 3). 



p = Me We -I- y^Wj-OT,- 

i 

4 = pip 

8. Results 



(112) 

(113) 
(114) 



We apply our ProDiMo model to a typical passive protoplane- 
tary disk of mass Mdisk - 0.01 Mg which extends from 0.5 AU to 
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1.5 2.0 2.5 3.0 3.5 
logTg[Kl 




r [AU] r [AU] 

Fig. 7. Gas temperature structure T^ir, z) in a model for a T Tauri type disk with Mjisk = 0.01 Mq (l.h.s.). See further model parameter in Table 5. 
On the rh.s. we show the results for the same parameter, if Tg = is assumed. The white dashed lines show Av = i and Ay = 10. 
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T„ = Ta assumed 
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Fig. 8. Dust temperature structure Ti(r,z)- The differences between the l.h.s. and the r.h.s. model are caused by the different density structures 
which are a consistent result of the entire coupled physical problem. 



500 AU. The central star is assumed to be a T Tauri-type "young 
sun" with parameters T^ff - 5770 K and L* = 1 Lq, and to emit 
excess UV of predominantly chromospheric origin as shown in 
Fig. 2. The stellar UV excess creates an unshielded UV radiation 
strength of about x - 2 x 10^ at 1 AU (see Eq. 41). Further pa- 
rameter of our model are summarized in Table 5. Our selection 
of elements and chemical species is outlined in Table 1, and the 
applied element abundances are listed in Table 3. 

The model uses a 150x150 grid of points which are arranged 
along radial and vertical rays which enables us to calculate the 
respective column densities and line optical depths in a simple 
way. The spatial resolution is much higher in the inner regions 
and the grid points are also somewhat concentrated toward the 
midplane. About half of the grid points are located inside of 
2.25 AU in this model to resolve the strong gradients in the ra- 
diation field and in the thermal and chemical structure occuring 
just inside of the inner rim. 



8.1. Disk structure 

The physical structure of the disk is a consistent result of all 
model components: dust radiative transfer, chemistry, and heat- 
ing and cooling balance. In order to explore how important the 
inclusion of the gas heating and cooling balance is for the re- 
sulting disk structure, we compare the full model (depicted on 
the l.h.s. of the following figures) to a comparison model (r.h.s.) 
where we have assumed T„ = throughout the disk. 



8.1.1. Thermal structure 

Figures 7 and 8 show the resulting gas and dust temperature 
structures of the models, respectively. The most obvious feature 
in Fig. 7 is a hot surface layer (Tg a; 4000 -7 000 K) which bends 
around the inner rim and continues radially to about 10 AU. This 
hot surface layer is situated above z/r ^ 0.13 in this model. Its 
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r [AU] r [AU] 

Fig. 9. Density structure n<H>('', z) as function of relative height z/r and log r. The red dashed line on the l.h.s. encircles hot regions Tg > 1000 K. 



aas in thermal balance 



r„ = T,i assumed 
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Fig. 10. Density structure of the puffed-up inner rim. Regions with ;j(h> < 10''^ are suppressed for output (shown in white color). The red dashed 
contour line on the l.h.s. encircles hot regions Tg > 1000 K. 



lower edge is not related to the vertical Ay but rather to the posi- 
tion of the shadow casted by the puffed-up inner rim. It coincides 
with the first occurrence of CO and other molecules like OH (see 
Fig. 12). The hot surface layer is optically thin, predominantly 
atomic (molecule-free) and directly heated by the stellar radia- 
tion in various ways (see Sect. 8.1.4). 

The shielded and cold regions in the midplane (z/r ^ 0.06) 
are characterized by small deviations between 7g and 7^, due 
to effective thermal accommodation between gas and dust. 
However, beyond some critical radius, here a; 100 AU, even the 
midplane regions become optically thin, and the interstellar UV 
irradiation causes an increase of Tg. We find midplane temper- 
atures up to Tg 27d around 400 AU in this model. The critical 
radius is related to Ay = 1 and increases with disk mass. 

The upper layers z/r ^ 0. 1 at r ^ 20 AU show no clear trend, 
both To < Ti and 7g > is possible, due to a complicated super- 
position of various heating and cooling processes. 

Apart from the thermally decoupled layers at the inner rim, 
the surface and the very extended layers, the disk temperature is 



mainly controlled by the dust continuum radiative transfer (see 
Fig. 8). 7d(r, z) shows all the features typical for protoplanetary 
disks (see e.g. Pascucci et al., 2004; Pinte et al., 2009). The mid- 
plane dust optical depth at 1 jjm is about 1.8 ■ 10^ in this model. 
The slightly different Td-results for the two models are caused by 
the different density structures (see Fig. 9) which depend on 7g. 
In case of the full model, the vertically extended inner regions 
scatter the star light and thereby heat the disk from above. 

8.1 .2. To flare or not to flare 

Figure 9 shows the resulting density structures of both models. 
The full model (l.h.s.) exhibits a remarkable vertical extension 
(up to z/r a: 1 ) of both the inner rim and the surface layers inward 
of 10 AU. According to Eq. (5), the vertical scale height H is 
approximately (assuming z<^r, cj- const) given by 



2rCj 



(115) 
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where H is defined as p(z) ~p(0) exp(-z^///^). The temperature 
ratio Tg/Td reaches values of about 10-30 in the hot inner rim 
and the surface regions, and since oc Tg, the disk is vertically 
more extended by about the same factor in comparison to the 
= Jd-model. This applies to the inner rim in particular, because 
it is hot even at z = 0, whereas the enhancing effect only starts 
at z/r ^0.1 in general. However, in the regions r ^ \ -7 AU, 
Tg is almost constant in the hot surface layer 4000 -7000 K) 
and so ///r increases further with increasing radius, and the disk 
reaches its maximum vertical extension here. 

It is noteworthy that the vertical density structure p(z) may be 
locally inverted. Since Eq. (5) is a pressure constraint, the den- 
sity must locally re-increase if Tg drops quickly with increas- 
ing height. This happens in the uppermost layers, in particular 
around 10 AU at z/r ^ 0.8, a region which causes the most nu- 
merical problems during the course of the global iterations. 

At larger radii ^ 30 AU, both models show a comparable 
vertical extension, characterized by a generally flaring structure. 
The "flaring" (increase of H/r with increasing r) is a natural con- 
sequence of the radial dust temperature profile varying roughly 
like Id(r) oc r with p ^ 0.25 in the midplane and p » 0.35 -0.45 
in the optically thin parts, so H/r oc r^. 



8.1 .3. The puffed-up inner rim 

Figure 10 shows a magnification of the density structure in the 
innermost regions. The figure demonstrates the large impact of 
the treatment of the gas temperature in the model on the result- 
ing disk structure. There is a rapid decline of the density be- 
tween n<H> = 10^ cm~^ and 10^ cm~^, which is caused by the 
steep Tg-increase at given pressure at the top of the shadow at 
z/r « 0.13 casted by the inner rim (l.h.s.). Therefore, such den- 
sities merely exist in the model close to the star, but the cool 
and dense midplane regions (> 10^ cm"^) are surrounded by an 
extended "halo" composed of thin hot atomic gas of almost con- 
stant density (n<H> = 10^ to 10^ cm~^) which extends as high up 
as z/r « 0.5. These results are astonishingly robust against vari- 
ation of the disk mass Mdisk between 10""* and 10"' M© — we 
always find the same kind of halo composed of the same kind of 
gas with the same densities. Only the midplane regions contain 
more or less cold matter, according to M^isk- 

The assumed position of the iimer rim at 0.5 AU in our model 
implies maximum dust temperatures of about 500 K, which is 
well below the dust sublimation temperature, and the shape of 
the inner rim is controlled by the radial force equilibrium at 
the inner edge which implies a smooth density gradient, see 
Sect. 3.1. In contrast, IseUa & Natta (2005) investigated the ef- 
fect of pressure-dependent sublimation of refractory grains on 
the shape of the inner rim. In reality, different kinds of refrac- 
tory grains will be present which have not only dififerent and 
pressure-dependent sublimation temperatures, but the dust tem- 
peratures are strongly dependent on dust kind due to dust opacity 
effects (see Woitke, 2006), which can be expected to result in a 
highly complex chemical structure of the inner rim. 

In comparison, the 7g = Tj-model does not possess the hot 
surface layers and, consequently, shows a much flatter structure. 
The inner rim is much less puffed-up causing the shadow border- 
line to be situated deeper. The iimer "soft edge" is likewise less 
extended, only from 0.5 - 0.61 AU in the Tg = 7d-model, whereas 
is extends from 0.5 -0.8 AU in the fuU model, or about 40% of 
the inner radius. 



8.1.4. Tliermal balance 

Figure 1 1 shows the most important heating processes (l.h.s.) 
and the most important cooling process (r.h.s.) in the full model 
of the disk with the gas being in thermal balance. Again, there is 
a clear dividing line at z/r -0.13 coinciding with the shadow of 
the inner rim, which separates the directly illuminated hot sur- 
face layers from the shielded and cold midplane regions. 

The central midplane of the disk below Ay 10 is dom- 
inated by thermal accommodation which assures Tg » (see 
also Kamp & Dullemond, 2004; Nomura & Millar, 2005; Gorti 
& HoUenbach, 2008). Since UV photons cannot penetrate into 
these layers, cosmic-ray ionization is the only remaining heat- 
ing process, mostly compensated for by thermal accommodation 
cooling. In the central midplane r ^ 1 AU, before H2O freezes 
out (see Fig. 12), there is additionally H2O rotational cooling, 
as well as some H2 quadrupole and CO rotational cooUng just 
below Ay 5= 10. 

Between Ay ~ 10 and z/r a; 0.13, the UV radiation can 
partly penetrate into the disk via scattering from above (see 
Fig. 4). This creates an active photon-dominated region with a 
rich molecular chemistry, where most of the abundant molecules 
Uke H2, CO, HCN, OH and H2O form, usually referred to as tiie 
"intermediate warm molecular layer" (Bergin et al., 2007). The 
layer is predominantly heated by H2 formation on grain surfaces 
and, with increasing height, by photo-efFect on PAH molecules. 
The gas temperature increases upward in this layer, e. g. from 
~ 200 K to ~ 700 K at 1 AU, but the additional heating can still 
be balanced by thermal accommodation in our model. 

The upper edge of the warm molecular layer is characterized 
by a thin zone of intensive CO ro-vibrational cooling. Above 
this zone, CO is photo-dissociated - below this zone, the CO 
lines become optically thick. It is this CO ro-vibrational cooling 
that can counterbalance the upwards increasing UV heating for 
a while, until the heating becomes too strong even for CO. This 
happens just at the upper end of the disk shadow z/r ^0.13 where 
the direct stellar irradiation becomes dominant. 

Above the CO layer, the temperature suddenly jumps to 
about 5000 K, all molecules are destroyed (thermally and radia- 
tively), and we enter the hot surface layer described in the pre- 
vious sections. This layer is predominantly heated by coUisional 
de-excitation of vibrationally excited H* (inner regions) and by 
PAH heating (outer regions). Although H2 is barely existent at 
these heights above the disk (concentration is 10"^ to 10"^, see 
Fig. 12), the few H2 molecules formed on grain surfaces can 
easily be excited by UV fluorescence, and these H2 particles un- 
dergo de-exciting collisions. This heating is balanced by vari- 
ous line cooling mechanisms. Since molecules are not available, 
atoms and ions like Oi and Fen are most effective. The non- 
LTE cooling by the wealth of fine-structure, semi-forbidden and 
permitted Fe I and Fe 11 lines has been investigated in detail by 
(Woitke & Sedhnayr, 1999), who found that in particular the 
semi-forbidden iron lines provide one of the most efficient cool- 
ing mechanisms for warm, predominantly atomic gases at den- 
sities n<H> = 10^ to lO''' cm"^ 

Since the stellar optical to IR radiation can excite most of the 
Fe n levels directly, radiative heating occurs. This "background 
heating by Fen" as referred to in Fig. 11 (l.h.s.) turns out to 
contribute significantly to the heating of the hot atomic layer 
close to the star (r ^ 10 AU). In fact, further analysis shows that 
the gas temperature in a large fraction of the hot atomic layer is 
regulated by Fpeii ~ Apeii, i. e. by radiative equilibrium of the gas 
with respect to the Fe n line opacity. Similarly, we find a small 
zone in the midplane just behind the inner rim where radiative 
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Fig. 11. Leading heating process (l.h.s.) and leading cooling process (r.h.s.) of the model in gas thermal balance. The black dashed contour line 
indicates an optical extinction of Ay = 10. 



equilibrium with respect to the water hne opacity is established. 
The regulation of the gas temperature via radiative equilibrium is 
a typical feature for dense gases in strong radiation fields, e. g. in 
stellar atmospheres. This behavior is rather unusual in PDR and 
interstellar cloud research from where most of the other heating 
and cooling processes have been adopted. 

The more distant regions ^ 50 AU are characterized by an 
equilibrium between interstellar UV heating (photo-effect on 
PAHs) and [Cii] 158 /zm, [Oi] 63 and 144 //m, CO rotational 
line cooling, and thermal accommodation (e. g. Kamp & van 
Zadelhoff,2001). 

8.2. Chemical structure 

The following discussion of the chemical results focuses on 
aspects that are relevant for an understanding of the two- 
dimensional disk structure. We restrict it to the most important 
atomic and molecular cooling species and the species that trace 
the dominant caiTiers of the abundant elements hydrogen, car- 
bon and oxygen throughout the disk (Fig. 12). A more detailed 
discussion of particular chemical aspects and their relevance to 
observations will be the topic of future work. 

8.2.1. Atomic and molecular hydrogen 

Inward of about 10 AU, the H/H2 transition occurs at the lower 
boundary of the hot surface layer. There is a very sharp gradient 
of UV field and gas temperature explained by the shadow casted 
by the dust in the inner rim. Above the shadow, the gas tempera- 
ture is high enough to efficiently destroy molecular hydrogen via 
H2 H- H — > 3H, and also by collisions with atomic oxygen. 



At larger distances, H2 can form on grain surfaces as soon 
as the dust temperature drops to about lOOK, where the forma- 
tion efficiency eiT^) increases sharply. This happens primarily in 
the secondary puffed-up regions around 10 AU. The formation 
of molecular hydrogen beyond this distance is mainly controlled 
by H2 self-shielding, which is an intrinsically self-amplifying 
(i. e. unstable) process. In addition, the gas density increases by 
a factor of ~ 2 when H2 forms at given pressure, which causes 
increased collisional H2 formation rates in comparison to the 
photo-dissociation rates. This H2 formation instability leads to 
local overdense H2-rich regions in an otherwise atomic gas at 
high altitudes at about 10 AU in our model. Other molecules like 
OH and H2O are also affected and these molecules can show 
even larger concentration contrasts as compared to H2 which 
causes the instability. 

8.2.2. Electron concentration and dead zone 

The electron density in the upper part of the disk is set by the 
balance between UV ionizations and electron recombinations of 
atoms and molecules. In the UV obscured, cold and icy midplane 
below z/r X 0.05, extending radially from just behind the in- 
ner rim to a distance of about 30 AU, the electron concentration 
drops to values below 10"*^, but cosmic ray ionizations main- 
tain a minimum electron concentration of ~ 10"'" throughout 
the disk, because the vertical hydrogen column densities in this 
model are insufficient to absorb the cosmic rays (A^H2 ~ lO^^cm"^ 
at r = 3 AU). An electron concentration of ~ 10"'*' is two orders 
of magnitude larger than the minimum value of ~ 10"'^ required 
to sustain turbulence generation by magneto-rotational instabil- 
ity (MRI), see (Sano & Stone, 2002). Thus, our model does not 
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Fig. 12. Ciiemical composition of tiie gas in a TTauri type protoplanetary disk with Mjisi^ = 0.01 Mq, showing the concentrations e, =n,/n(H). where 
iij is the particle density of kind ; and «(h> the total hydrogen nuclei density. Upper row: H, H2 and free electrons, second row: C^, C and CO, third 
row: O, OH and H2O and lower row: H2O ice, CO2 ice and CO ice. Note the different scaling for H, H2 and e" in the first row. 
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possess a "dead zone" in the planet forming region, which is 
different from studies about massive and compact, actively ac- 
creting disks (e.g. Ilgner & Nelson, 2006). 

8.2.3. C-", C, O, CO, OH, and H2O 

Outside the shadowed regions, the models clearly show the clas- 
sical C^/ C / CO / CO-ice transition as expected from PDR chem- 
istry (e.g. Kamp & Dullemond, 2004; Jonkheid et al., 2004; 
Gorti & Hollenbach, 2008). However, there are some important 
differences to note in the 1-10 AU range. The dominant form of 
carbon in the midplane is CH4. At those high densities, oxygen 
is locked up into HaO-ice, leaving carbon to form methane in- 
stead of CO. Above the icy regions, in a belt up to Ay « 10, 
water molecules evaporate from the ice and CO becomes again 
the dominant carbon and oxygen carrier. 

At radial distances between »0.8- 10 AU in the warm inter- 
mediate layer, the model shows a double layer with high con- 
centrations of neutral C, OH, H2O and other, partly organic 
molecules like CO2 HCN and H2CO (not depicted). This dou- 
ble layer is a result of the full 2D radiative transfer modelling in 
ProDiMo . The radial UV intensities drop quickly by orders of 
magnitude at the position of the inner rim shadow (z/r s; 0.13). 
The UV radiation field then stays about constant, until Ay ^ \ 
is reached, and also the vertical (-1- scattered) UV intensities de- 
crease. In combination with the downward decreasing gas tem- 
peratures and increasing gas densities, this produces two layers 
of hot and cold OH and H2O molecules with a maximum of C^ 
in between, van Zadelhofl" et al. (2003) have undertaken similar 
investigations showing that dust scattering leads to a a deeper 
penetration and redirection of the stellar UV into the vertical di- 
rection, with strong impact on the photo-chemistry. 

8.2.4. Ice formation 

The ice formation is mainly a function of kind, gas density 
and dust temperature. Hence, the location of the individual "ice 
lines" strongly depend on the disk dust properties assumed, such 
as total grain surface area, disk shape and dust opacity. Water and 
CO2 ice formation is mostly restricted to the midplane, where 
the densities are in excess of 10^*^ cm"^, the reason being mainly 
the reaction pathways leading to the formation of the gaseous 
molecules that form these ice species. In addition, UV desorp- 
tion counteracts the freeze-out of molecules in the upper layers 
at large distances from the star 

Inside 100 AU, densities are high enough to form water in the 
gas phase which subsequently freezes out onto the cold grains 
(Td £ 100 K). This is a consequence of our stationary chemistry 
that does not care about the intrinsically long timescale for ice 
formation (see Fig. 13). As densities drop and conditions for wa- 
ter formation in the gas phase become less favorable, oxygen 
predominantly forms CO, which freezes out at dust temperatures 
below ~ 25 K at large distances. There is an intermediate density 
and temperature regime (20-100 AU), where significant amounts 
of C02-ice are formed. 

8.3. Timescales 

An important question is whether our assumptions of gas energy 
balance and kinetic chemical equilibrium are vaUd in protoplan- 
etary disks. The cooling relaxation timescale is calculated as 



'''cool 
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where Q = r—A is the net heating rate. The chemical relaxation 
timescale is more complicated. We calculate it as 



Tchem = max Re{/1„} 

valid n ' 



(117) 



where A,, are the valid eigenvalues of the chemical Jacobian 
dFijdrij (see Eq. 64, without element conservation). Since F 
obeys A'ei auxiUary conditions (element conservation), there are 
A^ei redundant modes which have (mathematically) zero eigen- 
values. Numerics yields extremely small A for these modes, 
which must be disregarded. 

Figure 13 (l.h.s.) shows that the cooling timescale in the disk 
is smaller than typical evolution timescales by orders of magni- 
tude, and also smaller than typical mixing timescales (e.g. Ilgner 
et al., 2004), justifying our assumption of thermal balance. In 
particular, the gas is thermally tightly coupled to the dust via 
thermal accommodation in the midplane regions where the cool- 
ing timescale scales as Tcooi ~ 1 Ip- 

The chemical relaxation timescales (r.h.s. of Fig. 13) show 
that apart from the icy midplane, where Tchem can be as large as 
10^ yrs, the chemical relaxation timescale is typically 10"* yrs or 
shorter, and as short as ^ 1 - 100 yrs in the photon-dominated 
warm intermediate layer, where most spectral lines form. 

8.4. Gas emission lines 

In order to discuss from which part of the disk the various gas 
emission lines come from, we provide some simple estimates 
of column densities and excitation temperatures in this section. 
Full 2D non-LTE line transfer calculations will be covered in 
subsequent papers. 

The species column density A^thick required to achieve unit 
line optical depth can be calculated from Eq. (75). Assuming 
maximum and vanishing population in the lower and upper level, 
respectively (n; « ngp. 1m ~ 0), the result for r^" = 1 is 
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Table 6 shows some typical values of A^twck for permitted atomic 
resonance-lines, atomic fine-structure lines, and some rotational 
and ro-vibrational molecular lines. The actual critical column 
density is higher by a factor of Hsp/n; if the population of the 
lower level is less than maximum. 

Since the emission lines get saturated around t » 1, the ma- 
jority of the observable line flux originates from a surface region 
of thickness A^Jp' « A'thick- From the full model, we calculate the 
vertical species column densities A^™"^ and column density aver- 
aged gas temperatures (7'g)sp defineci as 



z 



(119) 



(rg>sp(r,iV-) = j^^J^^ J T,ir,z)ns,ir,z)dz (120) 



(116) 



Similar to Eq. (120), we define the column density averaged dust 
temperature (rd)sp. The results are shown in Fig. 14. Treating the 
column above A'^^p = A^thick as optically thin, ignoring deeper lay- 
ers for the line formation, and assuming LTE, the column density 
averaged gas temperatures {Tg)sp at depth N^^' = A^tMck provide 
an estimate of the expected excitation temperature of the observ- 
able fine flux. 
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Fig. 13. Cooling relaxation timescale Tcooi (l.h.s.) and chemical relaxation 



Table 6. Species masses in the disk and emission line characteristics. A 
velocity width of Ai'o = 1 km/s is assumed. 





mass [M(g] 


type 




Wih.ck [cm -] 


E„ [K] 


Cii 


0.055 


fine-struc. 


158 


2x 10" 


91 


Ci 


0.079 


fine-struc. 


370 


2x 10" 


62 


Oi 


4.6 


fine-struc. 


63 


3x 10" 


230 


CO 


4.2 


rotational 


400 - 2600 


-10'= 


-5-300 






ro-vib. fund. 


4.6 


-10" 


-3100 


H20 


0.007 


rotational 


50 - 300 


-10'^- 10'= 


-60- 1000 


H, 


2540 


ro-vib. 


2-28 


- 10^^ - 10^'' 


-500-8000 


Mgll 


0.15 


resonance k 


0.28 


6x 10" 


51300 



CII: The most simple case in Fig. 14 is the 157.7 //m fine- 
structure line of Cii. The column density reaches a maximum 
value of ~ 10'^ cm^, almost independent of radius r. This value 
falls just short of Mhick, meaning that the line remains mostly 
optically thin throughout the entire disk. The column of emit- 
ting C*" gas is situated well above the optically thick dust in the 
midplane as indicated by the missing Ay marks in Fig. 14. Since 
the outer regions have a much larger surface area as compared to 
the close regions, and the [Cii] 157.7 //m fine-structure line can 
easily be excited even in a cold low density gas out to 500 AU 
(£11 [K] = 91, «cr ~ 3 X lO"* cm"-'), the line flux is expected to be 
dominated by the outer regions. The [Cii] 157.7 fim line probes 
the upper flared surface layers of the outer disk. 

OI: The column of atomic oxygen gas responsible for the [Oi] 
63.2 yum fine-structure line extends deeper into the midplane 
and reaches column densities of ~ 10'^ cm"^, i. e. well beyond 
Mhick- Therefore, the line can be expected to be optically thick 
in most cases. The line saturates before a dust visual extinction 
of Ay = 0.1 is reached, i. e. also this line mainly probes the con- 
ditions above the optically thick midplane. Due to the higher en- 
ergy E„ required to excite the upper level, the line is expected to 
originate mainly from regions inward of about 100 AU. We esti- 
mate this radius roughly by the requirement that columns must 
provide temperatures 5 kTg ^ E,, (see Fig. 14, arrows on r.h.s.). At 
a column density of N^p - Nthick, the line intensity adopts an ex- 
citation temperature of ^ 40-70 K (r > 30 AU) and a; 500-1 000 K 
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timescale Tchem (r.h.s.). 



(r ^ 10 AU). In these inner regions, the [Oi] 63.2 //m line partly 
originates from the hot atomic surface layer of the disk (see 
Sect. 8.1.1) where {Tg)o peaks to about 1000-4000K at column 
densities ~ 10'*...10'^ cm"^, although eventually the excitation 
temperature lowers again as the column of oxygen gas extends 
into deeper, cooler layers. The final value of a: 500- 1000 K is 
~ 5 - 10 times higher than expected from Tg = models (see 
dashed lines in Fig. 14), stressing the necessity to include the 
gas energy balance in models for spectral interpretation. The 
[Oi] 63.2 fim line originates from the thermally decoupled sur- 
face layers inward of about 100 AU, above Ay ~0.l. 

CO: Both the lowest rotational line 1 ^ of CO as well as the 
CO fundamental ro-vibrational line (v, 7) = (1, 1) ^ (0, 0) need a 
column density of about 10'^ cm"^ to saturate. These CO column 
densities are exceeded by about 2.5 -4.5 orders of magnitude at 
all radii, i. e. these lines are optically thick throughout the entire 
disk. However, high-/ rotational or high-7 ro-vibrational tran- 
sitions test deeper layers because «/ <k «co and hence Mhick ^ 
10'^ cm"^. The observable line intensities are triggered mainly 
by the excitation criterion. For ro-vibrational transitions, for ex- 
ample, temperatures of the order of a lOOOK are only available 
inward of about 10 AU, where (Tg) has a strong negative gradi- 
ent, i. e. the emitting CO is situated just below the borderline be- 
tween the hot atomic and the warm intermediate disk layer, and 
we find values of about <rg)co(A^co = ^thick) ~ 1000 K. In the 
more distant columns, the ro-vibrational lines are also optically 
thick but the line source function is very small. The lower 
rotational lines until 5 — » 4 originate from the entire disk, but 
for higher rotational transitions, the emitting region shrinks re- 
markably. The CO rotational and ro- vibrational lines are mostly 
optically thick and probe very different regions in the disk, de- 
pending on the excitation energy E^. 

H2O: The situation for the rotational water lines is quite dif- 
ferent. Large amounts of H2O only form in deep layers, typi- 
cally below Ay - \. However, between 1 AU and 100 AU, the 
total H2O column densities are limited by about lO'^^-lO"" cm 
because of ice formation in even deeper layers. Since the rota- 
tional water lines have larger A„/ in general and get optically 
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Fig. 14. Column density averaged gas temperatures (r„> (Eq. 120. full lines) over vertical species column density for different radii labelled by the 
numbers in AU. The thin dashed lines show the column density averaged dust temperatures (T^i) for comparison. The small open circles, large open 
circles and large full circles mark vertical Ay = 0.1, Ay = I and Ay =10, respectively. The vertical arrows at the top of the figures mark the species 
column densities N^ck where the indicated lines become optically thick (Eq. 118). We add some horizontal (rightward) arrows for rotational and 
ro-vibrational lines to demonstrate that Aftuck is larger for high excitation lines. The leftward arrows on the r.h.s. of the plots indicate temperatures 
sufficient to thermally excite the upper level of the transitions 5 feT » £„. 



thick sooner, these column densities seem sufficient to saturate 
the low-lying rotational lines out to about 30 AU. Because gas 
and dust temperatures are coupled in the deep layers, the exci- 
tation temperatures of the rotational water lines in LTE with the 
local dust temperatures. Therefore, the Une-to-continuum ratio 
might be a problem for observations. 

Higher rotational lines are difficult to excite and the origin of 
such lines is probably limited to regions r^l AU. In particular, 
the inside of the inner rim is full of water and the temperatures 
here are high enough to excite a wealth of high-excitation ro- 
tational and probably also ro-vibrational lines, {Tg}H2o(N^'Q - 
A^thick) ~ 200 K- 2000 K between r = 0.6 AU and r = 0.8 AU. 



Another interesting region is the vertically extended zone around 
10 AU (see Figs. 9 and 12) which contains some amounts of hot 
water. We are currently investigating the impact of this hot water 
layer on the rotational lines as obversable with Herschel (Woitke 
et al., 2009). The rotational H2O lines probe the conditions in 
the midplane regions Av ~ 1 - 30, the inside of the inner rim, 
and possibly a vertically extended region with hot water around 
10 AU according to this model. 
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9. Conclusions and outlook 

This paper introduces a new code, ProDiMo, to model the phys- 
ical, chemical and thermal structure of protoplanetary disks. 
The strength of the new code lies in a fully coupled treatment 
of 2D dust continuum radiative transfer, gas phase and photo- 
chemistry, ice formation, heating & cooling balance, and the 
hydrostatic disk structure. In particular, we use the calculated 
radiation field as input for the photo-chemistry and as back- 
ground continuum for the non-LTE modelling of atoms, ions and 
molecules. The resulting gas temperatures determine the verti- 
cal disk extension, which in turn serves as input for the radiative 
transfer. Another advantage of the code is the robustness of its 
kinetic chemistry module which is applicable to densities be- 
tween 10^ to 10^^ cm"^; this makes possible to model complete 
disks ranging from about ~ 0.5 AU to 500 AU. 

Heating and cooling: The heating & cooling balance of the gas 
is the key to understand the vertical disk extension. The stellar 
UV irradiation, both direct and indirect via scattering, produces 
hot surface layers, in particular inside of ^ 10 AU, even without 
X-rays. We have included large non-LTE systems for Fe n and 
CO ro- vibrational fine transitions to counterbalance this heating 
in the disk regions with high p and high Tg, but find that these 
transitions also open new channels of radiative heating by the 
stellar optical -IR irradiation. The disk surface layers close to 
the star are in fact often stellar-atmosphere-like and character- 
ized by radiative equilibrium. More work is required to identify 
the important heating & cooling processes in these layers and 
the gas inside of the inner dust rim. 

Puffed-up inner rim and atomic lialo: Applying the new con- 
cept of "soft edges" to the inner rim of a T Tauri disk, the models 
show a highly pufFed-up inner rim extending up to z/r ~ 0.7, 
and an extended layer of hot (Jg * 5000 K) and thin (n(H> ~ 
10^-10^ cm~^) atomic gas reaching up to z/r =s 0.5 in the inner- 
most 10 AU. This "halo" is located above the intermediate warm 
molecular layer that surrounds the compact, dense, cold and icy 
midplane. It seems questionable that the highly pufFed-up struc- 
tures are hydrodynamical stable and further investigations must 
show how these features are related to mixing, winds and gas 
removal (Alexander, 2008; Woitke et al., 2008). 

Scattering and plioto-chemistry: The dust grains in the pufFed- 
up inner rim and the halo scatter the stellar UV light back onto 
the disk surface, which enhances the photo-chemistry and the 
photo-desorption of ice at larger distances. 

Chemistry: The surface regions of the model reveal the classi- 
cal PDR structure for H2, H, C+, C and CO. However, due to the 
full 2D UV radiation transfer in ProDiMo, a complicated multi- 
layered structure results for H2O and other organic molecules 
like CO2 and HCN, which depend sensitively on the model pa- 
rameters. The UV radiation field in the intermediate warm layer 
is reduced in two steps, first the pufFed-up inner rim blocks the 
direct path oF the radial photons From the star, and second the 
vertical and scattered photons are absorbed in the deeper layers. 
We find in particular two layers oFhot and cold water molecules. 
Mixing by hydrodynamical motions is likely to smooth out such 
structures (Ilgner et al., 2004; Semenov et al., 2006; Tschamuter 
& Gail, 2007). 

Cooling and chemical timescales: From the calculated relax- 
ation timescales we conclude that the assumption oF gas ther- 
mal balance and kinetic chemical equiUbrium should be suffi- 
cient For the interpretation oF most gas emission lines, although 
mixing may play a role. In the midplane, the chemical timescale 
is as long as 10** yrs due to ice Formation, but the spectral lines 



Form predominantly in higher layers where Tcooi ~ 1 ••• 100 yrs 
and Tchem ~ 1 ■■■ 10"* yrs, depending on distance. 
Emission lines: Disk emission lines originate mostly From the 
thermally decoupled surFace layers, where Tg > T^. A simple 
analysis oF the line characteristics and column densities shows 
that the [Cn] 157.7 //m fine-structure line probes the outer flared 
surFace oF the disk, whereas the [Oi] 63.2 jim line originates From 
slightly lower layers, but also probes the hot atomic gas inside oF 
10 AU. We identiFy this line as the key to test our models against 
observations. CO and H2O rotational lines probe the conditions 
in the intermediate warm molecular layer and. For low excita- 
tion lines, the outer regions. The origin oF most H2O lines is 
possibly restricted to regions ^ 30 AU, partly because the wa- 
ter abundances are higher there, and partly because the lines are 
thermally excited only in these layers. 

Observations: We intend to apply ProDiMo to a large sample 
oF observational disk emission line data that will be collected 
by the Pacs spectrometer on the Herschel satelUte, open time 
Key Program Gasps. Based on calculated chemical and thermal 
disk structures, detailed non-LTE line radiative transFer calcula- 
tions For the O i and C 11 fine-structure lines as well as some CO 
and H2O molecular lines will be carried out For analysis. We ex- 
pect to be able to determine the gas mass in disks From the line 
data, and to find tracers For hot inner layers. In the next decade, 
ProDiMo can be used to interpret Alma data which probe the 
physical conditions and the chemical composition oF the gas in 
the planet Forming regions oF protoplanetary disks. 
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